MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OP  STANDARDS -1963- A 


(I 


L  ?.&*'■ 


. 


r‘^*: 


00^' 


E& :; 


Nf’sss 


.ifej 


-‘•'“‘■ir’sS*, 

.<*98 


M  a  i  ■--'•■r4%.  , . 

pi  i  g  *«fsiiiigiiiil 


-S&jv 
i5r?  r<>S-f '.7:  V*’ 


..  .  ••  $ 


■*>*£  'l£’*‘”C  -V  ■  «’*  •- V*'  v«^"^W5;-','jf*^--«f  - ".'.  ,»►  .*.  ■>':«Ks:*iA<'‘'t"1’  *'V-  •«'-"-  J2.- •'■•' v>  g-  •'  -  -’•••-.  -A 

~  ■-  V  :  '  i  ^  -  ■:■ 

..r  :^;7:  J‘:p7- -’’-..gC-IV-"  -i  ■  '  ■ 


UNCLASSIFIED 


I  Lib1 


(tf  \  REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

ULBCaaMJuA**** _ _  ,  2.  GOVT  accession  no. 

RADcJrR-8^194  /  ^  \Ah-Ah  9f 

s.  RECIPIENT'S  CATALOG  NUMBER 

1, 91  . - 

j4.  title  fl&asAsnmy  " 

AN  IMPROVED  jr£IELD  SOLUTION  FOR  A  CONDUCTING 
flbDYOF  REVOLUTION-2^  /  P' 

*  * 

Whase  ^ep#t#u  S 

6-  performing  oyg.  report  number 

N/A 

J.  JktrTWORf.J  /  — 

Joseph  R./Mautz  j  [)K\ 

Roger  Fy4larrington  ^ 

*  F^62-79-C^|il , 

t  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Syracuse  University  - '  — ^ 

Department  of  Electrical  and  Computer  y 

Engineering,  Syracuse  NY  13210  ( fvj 

AREA  6  WORK  UNIT  NUMBERS 

11.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  v 

Rome  Air  Development  Center  (RBCT)  /It  ) 

Griff iss  AFB  NY  13441  1  ^  L 

14.  MONITORING  AGENCY  NAME  A  AOOfftSS(HNfffF*r*fil  Irom  Controlling  Olltco) 

Same 

Jun^lbsjr 

IS.  SECURITY  CLASS,  (ol  thia  report ) 

UNCLASSIFIED 

IS*.  OECL  ASSI  FI  CATION  '  DOWNGRADING 

n/Acheoule 

16.  DISTRIBUTION  STATEMENT  (of  thfg  Raport) 

Approved  for  public  release;  distribution  unlimited. 

IT.  DISTRIBUTION  STATEMENT  (a!  (A.  abetract  entered  in  Black  20,  II  different  tram  Depart) 

Same 

RADC  Project  Engineer:  Roy  F.  Stratton  (RBCT) 


»f.  KEY  WORDS  (Continue  on  rovoram  aid e  it  necaaaery  and  tdontlfy  by  block  number) 

Body  of  revolution 
Computer  program 
E-Field  solution 
Method  of  moments 

10  ABSTRACT  (Con limit  on  reeeree  elde  II  neceeearr  and  Identity  by  »loc*  number) 

The  electric  field  integro-diff erential  equation  for  electromag-? 
netic  scattering  from  a  perfectly  conducting  body  of  revolution  is 
solved  by  the  method  of  moments.  A  numerical  solution  is  obtained  by 
means  of  a  computer  program  which  is  described  and  listed.  This  com¬ 
puter  program  is  designed  to  handle  oblique  plane  wave  incidence 
efficiently.  Spatial  staggering  of  expansion  functions  for  the  ortho¬ 
gonal  components  of  the  induced  current  is  known  to  give  good  accuracy 


do  , ,  W3  COITION  OF  I  NOV  •*  IS  OBSOLETE 


_ UNCLASSIFIED 

SECURITY  CL  UNIFICATION  OF  TmI*  PAOE  ( 


(Pf*J 


(Cont'd) 

im  SwtErtwaR 


*H)l>  737 


T 


tern  20  (Cont'd) 

when  the  body  of  revolution  has  edges.  By  means  of  triangle  and  pulse 
expansion  functions,  this  spatial  staggering  is  achieved  without  the 
use  of  shifted  source  segments.  The  present  computer  program  is 
highly  competitive  with  other  available  programs  as  concerns  storage, 
execution  time,  and  accuracy. _ ^ 


K 


Ac  cession  For 
ST  IS  GftrtAI 
HOC  TAB 

'Ubamnooncad 

Justification - 


By  — - -  - 

|  PlstrifrnlA 

Availed'  1 ‘ ' 
jjUatl’ 

Wat.  ss't^al 


CONTENTS 


PART  ONE  -  SOLUTION  PROCEDURE  AND  NUMERICAL  RESULTS 


I.  INTRODUCTION - 


II.  METHOD  OF  MOMENTS  SOLUTION - 


III.  EVALUATION  OF  THE  MOMENT  MATRIX — 


IV.  EVALUATION  OF  THE  PLANE  WAVE  EXCITATION  VECTOR - 


V.  NUMERICAL  RESULTS- 


PART  TWO  -  COMPUTER  PROGRAM 


I.  INTRODUCTION - 


II.  THE  SUBROUTINE  ZMAT - 


III.  THE  FUNCTION  BLOG- 


IV.  THE  SUBROUTINE  PLANE - 


V.  THE  SUBROUTINES  DECOMP  AND  SOLVE- 


VI.  THE  MAIN  PROGRAM - 


REFERENCES - 


LIST  OF  TABLES 


Page 


Table  1.  Third  to  eleventh  arguments  of  ZMAT -  44 

Table  2.  Storage  of  matrix  elements  in  Z- -  50 

Table  3.  Fourth  to  ninth  arguments  of  PLANE -  58 


LIST  OF  FIGURES 

Page 


Fig. 

1. 

Triangle  function  T . (t) - 

5 

Fig. 

2. 

Pulse  function  P^ (t) - 

5 

Fig. 

3. 

Pulse  doublet  ~  T^ (t) - 

5 

Fig. 

4. 

Electric  current  for  axial  incidence  on  a  circular 
disk  of  radius  0.25A,  t  *  0  at  center - 

35 

Fig. 

5. 

Electric  current  for  axial  incidence  on  a  circular 
disk  of  radius  1.5A,  t  =  0  at  center - 

35 

Fig. 

6. 

Electric  current  for  axial  incidence  on  a  circular 
washer  of  inside  radius  0.4A  and  outside  radius  1.2A, 
t  =  0  at  inside  edge - - - 

36 

Fig. 

7. 

Electric  current  on  a  cone-sphere  of  cone  angle  20® 
and  sphere  radius  0.2A,  incidence  on  sphere - - - 

37 

Fig. 

8. 

Electric  current  on  a  cone-sphere  of  cone  angle  20® 
and  sphere  radius  0.2A,  incidence  on  tip - 

37 

Fig. 

9. 

Electric  current  on  an  open-ended  cylinder  of  radius 
X/(2tt)  and  length  X,  incidence  on  t  1  0 - — — - — 

% 

39 

Fig. 

10. 

Electric  current  on  a  spherical  shell  of  radius  0.2X 
with  axially  symmetric  aperture,  edge  at  t  *  0.471X, 
incidence  on  aperture - - - - - - — — - - 

39 

iv 


PART  ONE 


SOLUTION  PROCEDURE  AND  NUMERICAL  RESULTS 

I.  INTRODUCTION 

The  purpose  of  this  report  is  to  develop  an  efficient  numerical 
solution  to  the  E-field  integro-diff erential  equation  for  electromagnetic 
excitation  of  a  perfectly  conducting  body  of  revolution.  This  numerical 
solution  is  obtained  by  applying  the  method  of  moments  to  the  E-field 
equation.  The  E-field  equation  states  that  the  tangential  component  of 
the  total  electric  field  is  zero  on  the  surface  S  of  the  body  of  revolu¬ 
tion. 

The  problem  is  stated  in  Section  II  of  [1]  and  the  solution  is 
similar  to  that  in  Section  IV  of  [1].  Except  where  otherwise  indicated, 
the  notation  is  the  same  as  in  [1].  Equation  numbers  drawn  from  [1]  are 
preceded  by  1-.  For  instance,  (1-40)  denotes  equation  (40)  of  reference 
(11. 

The  following  differences  exist  between  the  present  solution  and 
that  in  [1].  In  the  present  solution,  the  approximation  to  the  generating 
curve  of  the  body  of  revolution  consists  of  half  as  many  straight  line 
segments  as  in  (1].  Otherwise,  the  t  directed  expansion  functions  are  the 
same  as  those  in  [1].  However,  for  <f>  directed  expansion-  functions,  the 
pulses  used  in  [2]  are  adopted.  Here,  t  is  the  arc  length  along  the 
generating  curve  and  <|>  is  the  azimuthal  angle.  The  testing  functions  are 
the  complex  conjugates  of  the  expansion  functions.  For  calculation  of  the 
elements  of  the  moment  matrix,  each  integral  with  respect  to  t'  over  each 
straight  line  segment  is  evaluated  by  using  n  -point  Gaussian  quadrature 


and  each  integral  with  respect  to  t  over  each  straight  line  segment  is 
approximated  by  sampling  at  the  midpoint  of  the  line  segment.  Although 
t  and  t’  are  both  arc  lengths  along  the  generating  curve,  t  denotes 
integration  over  a  testing  function  and  t'  denotes  integration  over  an 
expansion  function.  The  former  integration  is  called  a  field  integration, 
the  latter  a  source  integration.  As  in  [1],  n^-point  Gaussian  quadrature 
is  used  for  the  integration  with  respect  to  d> .  However,  the  method  [3]  of 
eliminating  the  singularity  is  used  to  fortify  the  Gaussian  quadrature 
integrations  with  respect  to  t'  and  4>  whenever  the  source  segment  is 
sufficiently  close  to  the  field  point.  For  calculation  of  the  elements  of 
the  excitation  vector,  n^-point  Gaussian  quadrature  is  used  for  the  t 
integration. 

With  regard  to  <p  directed  testing,  calculation  of  the  moment  matrix 
by  sampling  the  t  integrand  at  the  center  of  each  straight  line  segment  is 
equivalent  to  point  matching.  However,  for  t  directed  testing,  this  calcu¬ 
lation  can  not  be  viewed  as  simple  point  matching  because  each  t  directed 
testing  function  extends  over  two  intervals  and  therefore  must  be  repre¬ 
sented  by  two  Dirac  delta  functions  instead  of  one.  Furthermore,  the 
electric  charge  associated  with  each  t  directed  testing  function  is  also 
represented  by  two  Dirac  delta  functions. 

The  method  of  solution  formulated  in  Part  One  of  this  report  is 
Implemented  by  the  computer  program  described  and  listed  in  Part  Two.  The 
present  computer  program  takes  almost  twice  as  long  to  compile  as  that 
in  [1].  However,  for  axial  incidence  and  for  moment  matrices  of  roughly 
the  same  order,  the  present  program  with  nfc  ■  n^,  *  2  and  n^  *  20  executes 
almost  as  fast  as  that  in  [1]  with  N.  *  20.  For  moment  matrices  of  the 
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same  order,  the  present  computer  program  probably  executes  faster  than 
that  in  [2]  because  the  one  In  [2]  uses  twice  as  many  source  segments  and 
twice  as  many  field  points.  For  oblique  incidence,  several  moment  matrices 
are  required.  The  computer  program  in  [1]  calculates  the  moment  matrices 
one  by  one,  that  is,  each  moment  matrix  is  calculated  from  scratch.  How¬ 
ever,  the  present  computer  program  takes  advantage  of  the  fact  that  some 
intermediate  calculations  are  common  to  all  the  moment  matrices.  Hence, 
if  there  is  room  enough  to  store  all  the  moment  matrices  simultaneously, 
the  present  computer  program  should  execute  much  faster  for  oblique  inci¬ 
dence.  Results  obtained  from  the  present  computer  program  are  generally 
more  accurate  than  those  obtained  from  [1],  especially  for  bodies  of  revolu¬ 


tion  with  edges. 


II.  METHOD  OF  MOMENTS  SOLUTION 


The  boundary  condition  that  the  tangential  component  of  the  total 
electric  field  is  zero  on  S  is  expressed  by  (1-40)  and  supporting  equa¬ 
tions  (l-41)-(l-43)  .  Following  the  method  of  moments,  we  approximate  the 
electric  current  J  on  S -by 


l  (It.Jt.  +  I*jf.) 
nj-nj  nj-nj 


(1) 


and  substitute  this  J  into  (1-41).  In  (1),  JC .  and  J^,  are  known  expan- 

—  — nj  —nj 

sion  functions  and  I1",  and  1^.  are  unknown  coefficients  to  be  determined. 

nj  nj 

The  expansion  functions  JC.  and  J^.  are  defined  by 

-nj  -nj 


JC. 

-nj 


■Ht 


Jnj. 


j  =  1,2,  ...  P-2 
n  =  0 ,  +1 ,  +2 , . . . 


(2) 


J*  -  ».  ^ 

~nj  ~<P  Pj 


j  *  1,2,  ...  P-1 
n  =  0,  +1,  +2,... 


(3) 


where  and  u^  are  unit  vectors  in  the  t  and  <J>  directions,  respectively. 

The  j  which  appears  in  the  argument  of  the  exponential  in  (2)  and  (3)  is 

not  to  be  confused  with  the  j  which  appears  elsewhere  in  (2)  and  (3). 

The  former  j  is  /-T  and  the  latter  j  is  the  subscript  which  goes  from  1  to 

either  P-2  or  P-1.  The  function  Tj (t)  is  the  triangle  function  shown  in 

Fig.  1  and  p  is  the  distance  from  the  axis  of  the  body  of  revolution.  The 

function  P^(t)  *s  the  pulse  function  shown  in  Fig.  2  and  is  the  value 

of  p  at  t  *  t,  where  t.  is  the  center  point  of  the  domain  of  the  pulse, 
j  J 

The  purpose  of  the  scale  factor  1/p^  in  (3)  is  to  give  (3)  the  same  dimen¬ 
sion  as  (2),  namely,  1/length.  The  pulse  doublet  ^  T^ (t)  in  Fig.  3  is 


r 


used  later  on  In  the  method  of  moments  solution.  In  Figs.  1,  2,  and 
3,  t  is  the  arc  length  along  the  generating  curve.  It  is  assumed  that 
the  generating  curve  consists  of  P-1  straight  line  segments  where  P  is 
an  odd  integer  greater  than  or  equal  to  3.  The  jth  such  segment  extends 
from  t^  to  Its  length  is  .  The  expansion  functions  (2)  and  (3) 

are  especially  appropriate  if  the  body  of  revolution  is  an  infinitely 
thin  perfectly  conducting  surface  with  edges  at  both  ends  of  the  generat¬ 
ing  curve.  This  is  true  because  the  t  directed  electric  current  is  supposed 
to  approach  zero  at  an  edge  whereas  the  4>  directed  electric  current  might 
grow  large  there  [4], 

Testing  functions  4^  and  are  defined  by 


W* 

~~ni 


Ti(t)  -jn4> 
P  6 


i  =  1,2,  ...  P-2 
n  *  0,  +1,  +2,  . . . 


(4) 


«♦.  -  a*  —  e-j»* 

-ni  ^4»  p . 


i  =  1,2,  ...  P-1 
n  =  0,  +1 ,  +2 ,  ... 


(5) 


After  substitution  of  (1)  into  (1-41),  the  dot  product  of  (1-41)  is 
taken  with  each  testing  function.  These  dot  products  are  then  inte¬ 
grated  over  S.  As  can  be  derived  by  retracing  the  development  (1-40)- 
(1-65)  with  (1-46)  and  (1-47)  replaced  by  (2)-(5),  the  resulting  matrix 
equation  is 


> 

z* 

[t'l 

V! 

n 

n 

* 

n  I 

i 

Z* 

zn 

p 

V * 

n 

n 

1  n 
u  J 

n 

0,  +1,  +2,...  (6) 


where  the  Z  's  are  submatrices  and  the  I  's  and  V  ' s  are  column  vectors, 
n  n  n 
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The  matrix  of  the  Z  's  on  the  left-hand  side  of  (6)  is  a  square  matrix 
n 

called  the  moment  matrix.  The  column  vector  on  the  right-hand  side  of 

(6)  is  called  the  excitation  vector.  The  jth  element  of  I*  is  I*  and 

that  of  1^  is  1^  .  The  ith  elements  of  V12  and  are  given  by 
n  nj  n  n 

VtJ  =  -  [[  WC.  •  E^S  ,  i  =  1,2,  ...  P-2  (7] 

ni  n  JJ  ~ni  ~ 

S 

*  -  f[  •  ^dS  ,  i  -  1,2,  ...  P-1  (8! 

nr  n  J  J  ~ni 

S 

where  n  is  the  intrinsic  impedance  and  _E*  is  the  incident  electric 

field.  The  iith  elements  of  the  Z  's  are  given  by 
J  n 

*i+2  rj+2  2 

(z“)  .  =  j  I  dt  dt’  {k  T, (t)  T.(t’)(G  sin  v  sin  v’ 
n  ij  J  _  J_  1  3  5 

ci  ti 

+  G7cos  v  cos  v')  -  Gy  —  T^t)  Tj(t')}  (9 


.  .  fi+1  f  ,  d 

:n  }ij  =  ‘  J  dt  Pi(t)  j  dt'(k  p  Tj  (t’)GgSin  v'  +  nGy  ^  T.(t’)) 


fi+2  fj+l 


.  rltZ  N+l  ,  d 

)  -  —  I  dt  dt*  P  (t’Kkp'T^t)  G6sin  v  +  nG?  ^  T^t)) 

ij  Pj  J~  J~  J 


ci  ti 


**  -s  ri+1  rj+1  2 

(< \  -sirj  dtPi(t)  L  dt'pj(t')<k pp' 

1  3  t,  t. 


G5  -  n  G?) 


where 


°7  ’  G4  +  G5 


C  -jkR  „  , 

G4  =  2  d<}>  —  sin  (|-)  cos  (n<J>) 


d4>  — cos  <J>  cos  (n<j>) 


d$  — sin  $  sin  (n4>) 


R  *  /( p'-p)2  +  (z'-z)2  +  4pp’sin2  (|)  (17) 

Here,  k  is  the  propagation  constant,  p  is  the  distance  from  the  axis 
of  the  body  of  revolution,  z  is  the  rectangular  coordinate  along  this 
axis,  and  v  is  the  angle  that  the  tangent  to  the  generating  curve  makes 
with  the  z  axis.  The  angle  v  is  positive  if  p  increases  with  t  and 
negative  otherwise.  The  parameters  p,  z,  and  v  depend  on  t.  Their 
counterparts  p',  z',  and  v'  depend  on  t'.  The  ranges  of  values  of  i  and 
j  in  (9)-(12)  are  such  that  the  regions  of  integration  therein  move  from 
one  end  of  the  generating  curve  to  the  other  end.  It  is  understood  that 
n  -  °,  +1,  +2,  ...  in  (7)-(16). 

Note  that  the  quantity  defined  by  (14)  is  different  from  that 
defined  by  (1-62).  The  trigonometric  identity 

1=2  sin2(^)  +  cos  $  (18) 

was  used  to  express  (1-62)  as  the  sum  of  (14)  and  (15) .  Expression 
(14)  is  more  suitable  for  computation  than  (1-62)  because  the  inte¬ 
grand  in  (14)  is  always  finite. 
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III.  EVALUATION  OF  THE  MOMENT  MATRIX 


One  by  one  evaluation  of  the  elements  (9)-(12)  of  the  moment 


matrix  is  inefficient  because  of  the  overlapping  regions  of  inte¬ 


gration.  For  instance,  both  ^  and  contain  integrals 


tt. 


with  respect  to  t  over  the  ith  segment  (t.,  t.^,).  If  (Z  ).  ..  , 

i  i+1  n  l-l , j 


and 


.tt. 


(Z^  )  are  calculated  one  after  the  other,  these  integrals  must  either 


be  stored  or  calculated  twice. 

In  this  report,  the  contributions  to  (9)-(12)  are  accounted  for 
by  regions  of  integration  rather  than  by  matrix  elements.  Consider  the 
contributions  due  to  the  2-dimensional  region  of  integration 


t  <  t  <  t  . , 
p  -  -  p+1 


t  <  t’  <  t 
q  -  -  q+1 


This  region  of  integration  is  called  A^.  Integrations  in  (9)-(12)  are 


carried  out  over  A  for  {*  p}  or  possibly  {*  p  *}  ,  { ^  p  1 }  and  p  *}. 

pq  J=q  J=q  j=q-l  J=q-1 


For  all  other  values  of  i  and  j,  no  region  of  integration  in  (9)-(12) 

ri=P-1i  /i“P  1  f l-p-li  rimP) 


intersects  A  .  Setting  { .  p  ,}  ,  { .  p  . }  ,  { .  p  },  and  {  p}  successively 

pq  j«q-i  j=q-i  j-q  j*q 


in  (9)-(12)  and  counting  only  the  region  of  integration  A  ,  we  obtain 

pq 


fP+i  k+i 

e  j  dt  dt'{*  T.(t)T.(t')(G.sinvsinv'+  G.cosvcos  v')- 
n  ij  J  _  J  _  i  j  5  / 


(19) 


(Z^)  *  -  -i- 

"  *  : 


p+1 


dt  Pp(t)J  dt'(k  p  Tj  ( t ' )  G6sinv’+nG7  ^rT^t’)) 
Cq  (20) 


9 


The  asterisk  (  )  on  the  left-hand  sides  of  (19)-(21)  denotes  the  contri¬ 
bution  due  to  integration  over  the  region  A  .  Note  that  (22)  is  (12) 

pq 

with  ij  replaced  by  pq.  Because  (12)  has  no  overlapping  regions  of 
integration,  it  is  not  affected  by  the  change  from  calculation  by  matrix 
elements  to  calculation  by  regions  of  integration. 

Next,  each  integral  with  respect  to  t  in  (19)-(22)  is  evaluated  by 
using  the  approximation 

fP+1 

f(t)dt  -  f(t  )  A  (25) 

P  P 
t 

P 

where  f(t)  is  the  relevant  integrand  and,  as  indicated  in  Figs.  1,  2, 
and  3, 
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(26) 


Cp  =  2  (tp  +  tp+l) 


A  “  t  -  t 
P  P+1  P 


Application  of  (25)  to  each  integral  with  respect  to  t  in  (19)-(22) 
gives 


*  rq+I  2 

(ZtC),.  *  1A  dt'{k Y(t  )T.(t')(G_sin  v  sin  v'  +  G  cos  v  cos  v') 

n  ij  P  J_  ipj  5  P  '  p 


VV  j  p 


7  P 


'  VST 

P 


Z4*)  -  -  A  P  (t  )  fq  dt’ (k2T. (t')G,sin  v’  +  —  G?  T,(t'))  (29) 

n  pj  p  p  P  3  6  P  7  dt  3 


lZ\  -  A  P  dt'  P  (t')(~^—  T  (t  )G,sin  v  +  ~  G?  [-^  T  (t)]  ) 

n  iq  p  q  p  i  p  6  p  p  /  dt  i  t 


(Z**)  =  JA 

'  n  pq  p 


h+1  k2p*  n2 

P  (t  )  dt’P  (t,)(iVL-  G.  -  -S—  G7) 

P  d  n  D  5  DP  7 

v  -  M  o  a 


where  v  is  the  value  of  v  at  t  =  t  .  Incidentally,  v  =  v  for 
p  P  P 

t-  <  t  <  t-  because  the  generating  curve  was  assumed  to  be  straight 
P  P+1 

there.  In  (28)-(31),  Gj,  G6>  and  G?  are  given,  respectively,  by  (15), 
(16),  and  (13)  with  R  replaced  by  R  where 


in  (28)-(30)  is,  as  inherited  from  (19)-(21),  given  by  (23)  and  (24). 

Application  of  (25)  is  only  one  way  to  obtain  (28)-(31) .  Another 

way  to  obtain  (28) —(31)  is  by  approximating  the  G's  in  (19)-(22)  by 

their  values  at  t  *  t  .  This  amounts  to  immediate  rather  than  conse- 

P 

quential  replacement  of  R  by  Rp  in  (14)-(16).  A  third  way  to  obtain 
(28)— (31)  is  by  substituting  the  approximation 


Ti(t)  «  \  (Ai6(t-t.)  +  Ai+16(t-ti+1)) 

(33) 

Pp(t)  =  Ap6(t-tp) 

(34) 

&  Vc)  =  -  «<*-W 

(35) 

into  (19)-(22).  Here,  6(t)  is  the  Dirac  delta  function.  The  approxi¬ 
mation  (33)  preserves  the  value  of  the  surface  integral  of  the  t  com¬ 
ponent  of  the  t  directed  electric  current  (4)  on  the  portion  of  S  for 
which 

'p  - c  -  Vi  •  p " 1,J*  ••• P"1 

♦«  i  *  i  \ 

where  4)  and  <j>.  are  arbitrary.  Likewise,  the  approximations  (34)  and  (35) 
a  d 

do  not  alter  the  values  of  such  surface  Integrals  of  the  electric  current 
(5)  and  the  electric  charge  associated  with  either  (4)  or  (5). 

Equations  (28) -(31)  were  obtained  by  using  the  testing  functions  (4) 
and  (5)  and  invoking  either  the  approximation  (25)  or  the  set  of  approxi¬ 
mations  (33)-(35).  Can  a  set  of  effective  testing  functions  be  defined 
such  that  (28)-(31)  can  be  obtained  by  using  these  functions  and  no  auxili¬ 
ary  approximation?  Testing  functions  could  be  defined  by  substituting  (33) 


and  (34)  into  (4)  and  (5),  but  the  approximation  (35)  would  still  be 

required  in  order  to  obtain  (28) — (31) .  Unfortunately,  the  approximation 

(35)  is  not  consistent  with  the  approximation  (33).  Hence,  it  is  not 

possible  to  trace  (28) — (31)  to  effective  testing  functions. 

The  functions  P  (t1),  T.(t'),  777  T.(t'),  v',  and  p’  in  (28)-(31) 
q  j  at  j 

are  given  by 


P  (t’)  -  1 

q 


T  (t*)  -  f  + 


(-l)q“j(t'-t  ) 


t  (t»)  =  izll 

dt’  j*1  ;  A 


q-j 


,  j  -  q-l,  q 


»  j  -  q-i,  q 


v  *  v 


p’  =  p  +  (t’-t  )sin  v 

q  q  q 


(36) 

(37) 

(38) 

(39) 

(40) 


for  t  <  t'  <  t  .  Equations  (36)-(38)  can  be  obtained  from  Figs. 

q  q+i 

1,2,  and  3.  Equations  (39)  and  (40)  are  true  because  the  generating 
curve  is  straight  for  t  <  t'  <  t  Replacement  of  j,  q,  and  t'  by 

i,  p,  and  tp  in  (36) -(38)  gives 


p  (t  )  -  1 

(41) 

p  p 

VV  ■  7 

(42) 

[  **  T  ft)  ]  -  - — 

ldt  V  'Jt  A_ 

(43) 

Substitution  of  (36)-(43)  into  (28)-(31)  yields 
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fztc) 

(  n 


fq+i  k2 

jApJ_  dt'{T(1  + 

cq 


-)(Gcsin  v  sin  v  + 
5  p  q 


G.,  cos  v  cos  v  )  - 
7  P  q 


•*<ht  fq+1  k2  <-l)q_j2(t'-t  )  (-l)q'jnG7 

(Z9C)  -  -  A  dt'&-  a  +  - 7 - Mg, sin  v  +  ——7 - -) 

n  pi  p  2  A  6  q  pA 

q  P  q 

q  (45) 


fq+l  Jt  (f-t  ) 


fctd>  rqri 

r\n  -  An  df  (-V  a  + 

n  iq  P  J_  2 


(-1)H  "nG, 


- sin  v  )G,sin  v  +  - r - ) 

P  q  6  p  p  A 

q  M  r  qp 


(z  =  jA 

n  pq  p 


(q+1  2 
dt'(k  (1  + 


(t’-t  )  2 

— - — sin  v  )G  -  -~n  —  G  ) 
p  q  5  p  p  7 

q  M  P  q 


Equations  (42)- (45)  are  rewritten  as 


n  j  k  A  A 

(Z  C). .  =  - (G_  sin  v  sin  v  +  G,  cos  v  cos  v  )  + 

n  i]  8  '5a  p  q  7a  p  q 


(-l)q*3jk2A  A 


(G_.  sin  v  sin  v  +  G_,  cos  v  cos  v  )  - 
5b  p  q  7b  p  q 


i  s. 


i 


2  2 
k  A  A  sin  v  .  k  A  A  sin  v  nA^ 

<<  >pj  ’  -  ^ - 3>°6a  -  1 H - 1>=6b  *  (^S.) 

FJ 


<%„  ■  < 


k  A  A  sin  v 


A  sin  v  .  nA 

1)(G6a  + G6b>  +  V>G7a 


2 

. .  k  A  A  A  sin  v  nA„  nA^ 

«*>„  -  2i  h-Mho*.  *  — 2  o,k>  -  (^X^c,  ) 


n  pq 


2p  5b'  2p  2p  7a 
q  <1  P 


where  i  is  either  p-1  or  p  and  j  is  either  q-1  or  q  and  where 


2  fq^ 

G  =  G  dt' 

ma  A  m 


in  *  5,6,7 


2  ^q+l 


2  L  fq 

-<r> 

n  *  — 


t  )G  dt’ 
q  m 


Equation  (13)  is  used  to  rewrite  (52)  and  (53)  as 


G_  -  G.  +  G. 
7a  4a  5a 


G7b  *  G4b  +  G5b 


c 

2  f*1 

O  ^  — 


G  (t ’  -  t  )dt ' 
m  q 


m  *  4,5,6 


2  Wl 


<f)  r 

o  J  — 


(f  -  t  )G(t» 
q  m 


t  )dt' 

q 


The  argument  (t'  -  t^)  supplied  with  Gm  in  (56)  and  (57)  comes  into  play 

later  on.  Substitution  of  R  for  R  in  (14)-(16)  produces 

P 

it  -jkR 
(  P  »  , 

G4(t’  -  tq)  -  2  I  d<j>  -  sin  (|o  eos(n<|))  (I 

0  P 


it  -jkR 

(  e  p 

Gj(t 1  -  t  )  =  d(J>  -  cos  <j)  cos(n<}>)  (59) 

0  P 

ir  -jkR 

f  e  p 

Gg(t '  -  t  )  *  d4>  -  sin  4>  sin(n4>)  (60) 

0  P 

where  is  given  by  (32).  In  (32),  p'  is  given  by  (40)  and  z*  by 

z’  =  z  +  (t*  -  t  )cos  v  (61) 

q  q  q 

Equation  (61)  is  true  because  the  portion  of  the  generating  curve  for 
t  <  t'  <  t  ,,  is  straight. 

q  q+i 

Evaluation  of  the  integrals  in  (56)  and  (57)  by  means  of  an 
nt~point  Gaussian  quadrature  formula  gives 


n 

t 

(n  ) 

i  <■>,) 

II 

oe 

A£' 

G.(2\*V  > 

m  *  4,5,6 

(62) 

nt 

(n) 

(n  )  .  (n  ) 

G  .  -  l 
mb  vLml 

A£’ 

*£’  Gm(lV£’  > 

(63) 

(nc)  (nt) 

where  the  abscissas  x^,  and  weights  A^,  are  tabulated  In  Appendix  A 

of  [5]  for  several  values  of  nt-  Application  of  an  n^-point  Gaussian 

quadrature  formula  to  the  integrals  in  (58)-(60)  and  replacement  of 
1  (nt) 

(t'  -  t  )  by  r  A  x. ,  result  in 
q  z  q  x. 


1  <«,>  "♦  K>  e"jkRpr*  2  h 

GAj  Ao  V  >  ■  v  l  kl  $  "Tr -  sin  (*2  )c°8 (**♦#)  (64) 

4  2  q  X,  1-1  *  *%£’*  L  % 

,  (n. )  n4>  (n.)  "jkRPri 

G5(2  Aq  X£’  >  "  2  l  A£  kR  —  '  C0S  V08(n*£>  (65) 

£-1  p £'£ 
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> 

i 


where 


where 


.  (n  )  4> 

G6(2  Aq  X£’  >  =  f  I  A: 


<y 


kR 


pi' l 


sin  ((i^sinCntfi^) 


P* 


•i  ■  + 


2  ^2, 

4p  p'sin  (-=-) 
P  P  2 


1  (nt) 

P’  =  P  +-A  x„ ,  sin  v 
q  2  q  V  q 

1  (nt) 

z '  =  z  +  —  A  x  „ ,  cos  v 

q  2  q  V  q 


(66) 


(67) 


(68) 


(69) 


h  =  \  (x2  *  +  2) 


(O 


(70) 


1  t 

Calculated  values  of  G  (-  A  x. ,  )  from  (64)-(66)  are  substituted 

m  2  q  Jc 

into  (62)  and  (63)  in  order  to  evaluate  G  and  G  ..  The  resulting 

ma  mb 

values  of  G  and  G  ,  are  then  substituted,  either  directly  or  through 
m3  mb 

the  intermediary  equations  (54)  and  (55),  into  formulas  (48)~(51)  for 
the  elements  of  the  moment  matrix. 

The  values  nf  =  2  and  n^  =  20  are  suggested  whenever  the  field 
point  is  not  close  to  the  source  segment.  If  the  field  point  is  close 
to  the  source  segment,  the  method  of  eliminating  the  singularity  [3] 
is  used.  Since  double  integrals  are  involved,  three  variations  of  the 
method  are  possible.  These  variations  are  called  methods  1,  2,  and  3. 

In  method  1,  elimination  of  the  singularity  is  applied  to  the  integra¬ 
tion  with  respect  to  t'.  In  method  2,  elimination  of  the  singularity 
is  applied  to  the  integration  with  respect  to  <p.  In  method  3,  it  is 
applied  to  the  double  integral.  In  methods  1  and  2,  the  "singular  part" 
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of  the  integrand  is  subtracted  out,  numerical  integration  of  the 
resulting  finite  integrand  is  performed  with  respect  to  one  of  the 
variables,  the  integral  (with  respect  to  this  variable)  of  the  "singular 
part"  is  added,  and  then  numerical  integration  with  respect  to  the  other 
variable  is  done.  In  method  3,  the  singular  part  of  the  integrand  is 
subtracted  out,  numerical  integration  of  the  resulting  finite  integrand 
is  performed  with  respect  to  both  variables,  and  then  the  double  integral 
of  the  "singular  part"  is  added.  Method  3  is  preferable  to  either  of 
methods  1  and  2  because  the  final  numerical  integration  in  methods  1  and  2 
may  involve  a  singular  integrand.  However,  if  what  is  deemed  to  be  the 
"singular  part"  can  be  integrated  analytically  with  respect  to  only  one 
of  the  variables,  then  either  method  1  or  method  2  is  applicable,  but 
method  3  is  not. 

Use  of  method  1  is  now  demonstrated.  From  (56)-(60),  the  required 
integrals  with  respect  to  t'  are 


G 

a 


dt' 


t 


q 


dt' 


(71) 


(72) 


The  above  expressions  are  rewritten  as 


where 


al 


ff 

a  J  — 


t  -jkR 
q+1  e  P-1 


kR 


dt’ 


(75) 


2  fq+1  dt’ 


a2  A 


kR 


(76) 


2  2  P+1  e"jkRp-l  , 

Gbi  -  (t’  -  V(n5r-: ">dt 

q  p 

q 


(77) 


Gb2  - 

q 


2  rq+1  Ct  *  -  t  )dt' 


kR 


(78) 


Application  of  n^point  Gaussian  quadrature  to  the  right-hand  sides  of 
(75)  and  (77)  gives 

n 


t  (n  ) 

G.i '  *4  v G 


(79) 


t  (nt)  (nt) 

5bl  =  x£’  V  G 


(80) 


where 


-jkR. 


kR 


kR 


kR 


p-;L  -  sin  {.~2>  (sin(-jt)  +  j  cos  (-^v) 


kR 


kR 


(  2  } 


(O 


(81) 


1  t 

where  R  is  to  be  evaluated  at  (t*  -  t  )  ■  -r  A  x.,  .  The  purpose  of 

P  q  2  q  X 

the  alternate  form  of  G  on  the  extreme  right-hand  side  of  (81)  is  to 
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To  reduce  roundoff  error,  (88)  is  rewritten  as 


where  G ■  is  to  be  evaluated  at  <|>  =  <p^  given  by  (70).  Equations  (91)- 

(93)  are  also  valid  with  a  replaced  by  b.  Calculation  of  and  G^ 

should  be  according  to  the  development  (73)-(90)  only  for  those  values 

of  4>„  for  which  r  is  either  smaller  than  or  comparable  to  A  .  If  r 
i  pq  q  pq 

is  considerably  larger  than  A^,  pure  Gaussian  quadrature  is  adequate. 

Use  of  method  2  is  now  demonstrated.  Since  the  integrands  of 
(58)  and  (60)  are  fairly  well-behaved,  method  2  is  applied  only  to  (59). 
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In  method  2,  G,  and  G  are  calculated  according  to  (62)  and  (63)  with 

33  Ju 

1  (nt> 

Gt(r  A  x„,  )  given  not  by  (65)  but  by 

5  2  q  x. 


,  (n  )  _  ^  (n  )  jkV’!L 

S(2  }  "  2  A£  “  kRM7A  '  cos  ^4  cos(n<t,j^)  - 


“I  l 

^  n  _ 


k  /(p’-pp)2  +  (z'-zp)2  +  p'pp4>2 


P=0  k  ^/(p'”Pp)2  +  (z’-zp)2  +  p'p^ 


From  formula  200.01.  of  Dwight  [6], 


■0  k  / (p'-p  )2+(z'-z  )2  +  p'p  ?  k^p 


- -  log  (u  +  Vl+u2) 


where 


ir/p'p 

_ E _ 

f(p'-P„)2  +  (*’-*  )2 

P  P 


Equation  (94)  should  be  used  only  for  those  values  of  t'  for  which  p^ 

I  2  2 

is  considerably  larger  than  y(p'-pp)  +  (z'-zp)  .  Otherwise,  the  pure 
Gaussian  quadrature  of  (65)  is  adequate. 

Use  of  method  3  is  now  demonstrated  for  the  case  in  which  p*q. 
Method  3  is  applied  only  to  the  calculation  of  G^a  because  G^a  is  the 
only  integral  in  (56)  and  (57)  whose  integrand  is  not  bounded.  We  write 
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7T  ♦  <“*>  11  (nt)  e  p* 

G.  *  7  l  A.  9  cos  4>.cos(n(J>  )  £  A  ,  -rg - 

53  ZW  *  *  £'=1  p£'£ 


n<j)  (n  )  nt 

- 1  l  *j  *  I 


£' 


£'=1  /  A  (n  )  2 


2  ( 

+  f  d* 

n  ^  ^  - 


kV('fx£It)  +  V* 


2.  2,2 


2.2  4  0  t  k/(t’-t  )  +p  4> 


q  q 


Because  of  the  formula 


[x  log  (y  +  Jx2+y2)  +  y  log  (x  +  Jx2+y2)  ] 


the  double  Integral  in  (97)  is  tractable. 


f  f  d»  f+1  ,  -t'2  2-2  ■  vtr  ><■«  [: 

q  o  t“  kj(t'-t  >2  +  pV  qL 


A  +  V  1  +  <  A  >  j  + 
q  q  J 


2irp  r  A  /  A  _ 

+  ^l042-^Wl+O2 

q  L  q  q 


In  e3ch  of  methods  1,  2,  and  3,  an  attempt  is  made  to  subtract 
out  the  singularity  due  to  1/R^  in  (58)-(60).  In  method  1,  1/Rp  itself 
is  subtracted  out.  In  method  2,  the  approximation 


p'-p  )2  +  (z’-z  )2  +  P_P'4>2 
V  p  P  P 

to  1/R  is  subtracted  out.  For  comparison,  R  is  given  by  (32).  In  method  3 
P  P 


the  approximation 


f(p'-Pp)2  +  (z’-Zp)2  +  PpPq4>2 


—  UI»".MWP>IF 


to  1/Rp  is  subtracted  out  for  p=q.  Because  the  double  integral  of 
this  approximation  is  tractable,  method  3  can  be  extended  to  cover 
the  case  in  which  p  f  q.  For  p  t  q,  the  alternate  approximation 

1 


n/(p'-p  )2  +  (z'-z)2  +  P_P, 

*  p  p  pi 


p  min 


to  1/R  merits  consideration.  Here,  p  .  is  the  value  of  p’  at  that 
p  min 

2  2 

value  of  t'  which  minimizes  (p’-p  )  +  (z'-z  )  .  No  matter  which  of 

P  P 

the  above  two  approximations  to  1/R^  is  used,  the  closed  form  expres¬ 
sion  for  its  double  integral  is  rather  complicated  and  vulnerable  to 
roundoff  error.  For  this  reason,  method  3  was  used  only  for  p=q. 

For  p  ^  q,  the  decision  whether  to  use  methods  1  or  2  is  based 

on  comparisons  of  A  with  d  and  p  with  d  where  d  is  the  distance 
q  o  q  o  o 

from  the  field  point  at  t  =  tp  to  the  nearest  point  on  the  qth  source 
segment.  The  distance  between  the  field  point  at  t  =  t^  and  the  point 
(t’,<}>)  on  the  qth  source  segment  is  given  by  (82)  or  (83).  It  is 
evident  that  the  minimum  of  (82)  occurs  at  <t>  =  0  because  neither  p^  nor 
p'  of  (40)  can  be  negative.  At  $  ■  0,  (84)  and  (85)  specialize  to 

t  =  (p  -  p  )  sin  v  +  (z  -  z  )  cos  v  (100) 

o  q  p  q  q  p  q 

d*  =  I (p  -  P  )  cos  v  -  (z  -  z  )  sin  v  I  (101) 

1  q  p  q  q  p  q 

The  asterisk  (  )  on  the  left-hand  sides  of  (100)  and  (101)  indicates  < 

that  <}>  ■  0.  Minimizing  (83)  with  respect  to  t'  on  the  qth  source 

A  Aq 

segment  where  — ^  <  t1  -  t  <  -jS  we  obtain 
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(102) 


~N 

p  j4  q 


A  <  d 

t  q  —  o 


> 


c 


,P  <  d  i 
<fe  q  —  o  \ 


Case  1 

Pure  quadrature 


(103) 


then  the  pure  quadrature  of  (62)  -  (66)  is  used  to  calculate  G  and 


G  ,  .  Here,  c  and  c,  are  constants  for  which  the  values 
mb  t  <J) 


Ct  = 


are  suggested .  If 


c ,  p 

<p  q 


<  d  \ 
—  o  ' 


then  method  1  is  used .  If 


0.1 


P  +  q 

4  c  A  >  d 
2  t  q  o 


Case  2 
Method  1 


P  J*  q 

c ,  p  >  d 

<P  q  o 


p  =  q 


Case  3 
Method  2 


then  method  2  is  used.  If 

Case  4 

Methods  1  and  3 
then  both  methods  1  and  3  are  used. 


(104) 


(105) 


(106) 


(107) 
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The  strategy  in  (103)  and  (105)-(107)  is  based  on  the  assumptions 


l 


that  the  Gaussian  quadrature  integration  with  respect  tot' must  be 
fortified  only  when  A^  is  large,  and  that  the  Gaussian  quadrature  inte¬ 
gration  with  respect  to  4>  must  be  fortified  only  when  is  large.  The 
integration  with  respect  to  t'  could  not  be'  fortified  for  ct^q  >  *-n 

Case  -3  because  methods  1  and  2  can  not  be  applied  simultaneously  and 


because  it  was  decided  earlier  to  limit  use  of  method  3  to  Case  4.  How¬ 
ever,  pure  Gaussian  quadrature  should  still  give  a  fairly  accurate  evalu¬ 
ation  of  this  integral  with  respect  to  t'  because  of  the  following  reason¬ 
ing.  Since  p'  is  large,  difficulty  can  only  occur  when  0  is  small.  Further- 


It  is  evident  that  A  <  2d  if  p  j*  q,  if  all  the  A  are  equal,  and  if  the 

q  —  o  q 

generating  curve  does  not  fold  back  on  itself. 
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IV.  EVALUATION  OF  THE  PLANE  WAVE  EXCITATION  VECTOR 

Consider  the  elements  (7)  and  (8)  of  the  excitation  vector  for  a 
9-polarized  incident  plane  wave  defined  by 


i  t  "jVX 

E  =  Ugkn  e 

and  also  for  a  ^-polarized  incident  plane  wave  defined  by 

.  .  -jk  •  r 

i  t  J — t  — 

E  =  u'Tkri  e 


(108) 


In  (108)  and  (109), 


k  =  -  k(u  sin  0  +  u  cos  0  ) 

— t  —x  t  ~z  t 


Ufl  =  u  cos  0  -  u  sin  0 

—0  —x  t  — z  t 


(109) 


(110) 


(111) 


u*  =  u  (112) 

~y 

where  0^  is  the  angle  of  incidence  and  where  u  ,  u  ,  and  u  are  unit 
vectors  in  the  x,y,  and  z  directions,  respectively.  Also,  jr  is  the 
radius  vector  from  the  origin.  The  origin  must  lie  on  the  axis  of  the 
body  of  revolution  because  this  axis  is  the  z  axis.  Substitution  of 
(4),  (5),  and  (108)  into  (7)  and  (8)  gives 


tft  n 

V  *  j  7rk  dt  T.  (t){j  sin  v  cos  0  (J  --J  .)  -  2  cos  v  sin  0  J  }e 

ni  J  l  t  n+i  n— i  t  n 


jkz  cos  0 


Vfl  '  ]"’k  J.  dt  t  VC>  <WJn-l>“8  9t 


jkz  cos  0 


(113) 


(114) 


where  VC®  is  VC.  for  E*  given  by  (108)  and  V^®  is  V^.  for  £*  given  by 
ni  ni  m  m 

(108).  In  (113)  and  (114), 


J  =  J  (kp  sin  0.) 
n  n  t 


(115) 


where  is  the  Bessel  function  of  the  first  kind.  Likewise,  substitu¬ 
tion  of  (4),  (5),  and  (109)  into  (7)  and  (8)  gives 

.  ri+2  jkz  cos  0 

ni  *  "  j  7rk  J  dt  Ti(t)(Jn+l'+  Jn-l>Sin  V  6  (116) 


V** 

ni 


1"+1'k  L  dt  t  pi(t)  <Vi 


n-l 


jk2  cos  0 
)e  1 


(117) 


where  the  second  superscript  on  V  on  the  left-hand  sides  of  (116)  and 
(117)  denotes  excitation  by  the  <j)-polarized  incident  plane  wave  (109). 

The  manipulations  required  to  obtain  (113)-(117)  are  similar  to  those 
used  in  the  derivation  of  (1-95) . 

The  contributions  to  (113)  and  (116)  due  to  integration  with 
respect  to  t  from  tp  to  tp+1  are  expressed  by 

*t0  n  f  P*1  jkz cos  0 

V  .  ■  j  irk  dt  T  (t){  j  sin  v  cos  8  (J  ,--J  ,)~2  cos  v  sin  0kJ  }e 

ni  J_  i  t  n+i  n-l  t  n 

CP  (118) 


*t<b  n  W  jkz  cos  0t 

Vni  "  3  nk  J  dt  Ti(t)(Jn+l+Jn-l)sin  V  6  <119) 

t” 

P 

A 

where  i  is  either  p-1  or  p.  The  asterisk  (  )  on  the  left-hand  sides  of 
(118)  and  (119)  denotes  the  contribution  due  to  integration  from  tp  to 
t . . , .  First,  v  is  replaced  by  v  in  (118)  and  (119).  Throughout  (114), 

p+i  p 

(117),  (118),  and  (119),  P^t),  T^(t),  and  p  are  expressed  according  to 
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(36),  (37),  and  (40),  respectively.  Then,  i  is  replaced  by  p  in  (114) 
and  (117)  to  make  those  equations  compatible  with  (118)  and  (119).  The 
results  of  the  above  substitutions  are 


*t0  inTrk  f^1 
Vni  "  V"  L  dt  (1  + 


(-l)P_i2(t-t  ) 


A - *-Hi  sin  v  cos  VVl'-W  - 

P 

jkz  cos  0 

-2  cos  v  sin  0  J_)  e  (120) 

P  t  w 


V*e  -  Ak 
np 


fP+1 


(t-t  )sin  v  jkz  cos  0 

dt  (1  ♦  — E- - E><-V1«,-l)c°s  0t  ' 


(121) 


*td»  .  jnTTk  fP+1 
ni  "  *  2 


dt(l  + 


(-l)p_i2(t-t  ) 


)(J  +J  )sin  v  e 
n+1  n-1  p 


jkz  cos  0 


(122) 


,,({><|>  .n+1  i 

VYY  *  j  irk 
np 


P+1 


(t-t  )sin  v 


dt 


(1  t  - t- - E)<Jn+1-Jn-l>  ' 


jkz  cos  0 


(123) 


where  i  is  either  p-1  or  p  in  (120)  and  (122) . 

Equations  (120)-(123)  can  be  rewritten  as 


.  a  j1^^  irkA  sin  v  cos  0  jniTkA  cos  v  sin  0 

$t0  »  - - E - E - 1  (F  _F  J - e_ - E - 1  p  + 

4  n+1, a  n-l,a;  2 


ni 


na 


jn+^irkA  sin  v  cos  0_  jnirkA  cos  v  sin  0 

+f-llp~i{- _ 2 _ 2 _ -  (F  -F  ) _ 2 - 2 - -  F  > 

1  '■  '  n+l,b  n-l,b' 


nb 
(124) 


A  sin  v 


,a  j  TTkA  COS  0  -  _ 

v*e - * - -c  «pn+1>  «n.1>a)  ♦ 


np 


(125) 


1 

1 
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where 


h 


sin  v 

P 


COS  V 

p 


(134) 

(135) 


The  calculation  of  the  plane  wave  excitation  vector  would  be  most  nearly 
consistent  with  the  calculation  of  the  moment  matrix  if  n^,  =  1.  For 
nT  =  1,  the  extra  data  x^^  =  0  and  =  2  must  be  supplied.  Now, 

assuming  that  n  >  1,  it  could  be  said  that  nt~point  quadrature  is  more 
accurate  than  1-point  quadrature.  The  n£ -point  quadrature  data  are 
already  available  because  they  were  used  to  calculate  tne  elements  of 
the  moment  matrix  in  Section  III.  With  nt  fixed  at  2,  results  were 
calculated  for  both  *  1  and  n^  =  2.  It  was  difficult  to  tell  which 
results  were  more  accurate.  The  numerical  results  presented  in  Section  V 
were  obtained  by  using  n^  *  n^  =  2. 


V.  NUMERICAL  RESULTS 


Computer  program  subroutines  have  been  written  to  calculate  the 
elements  of  the  moment  matrix  and  the  elements  of  the  plane  wave  exci¬ 
tation  vector.  These  subroutines  are  described  and  listed  in  Part  Two  of 
this  report.  They  were  used  to  calculate  the  electric  currents  induced 
by  a  plane  wave  axially  incident  on  two  circular  disks,  a  thin  washer, 
a  cone-sphere,  an  open  cylinder,  and  a  spherical  shell  with  an  axially 
symmetric  aperture.  The  magnitudes  of  these  electric  currents  are 
plotted  in  this  section. 

For  axial  incidence,  0  is  either  0°  or  180°  and  the  only  non¬ 
zero  excitation  vectors  for  the  0-polarized  plane  wave  (108)  are 


—  — 

-  n 

-1 

s 

V1 

i 

V* 

?> 

t 

-1 

L 

It  is  evident  from  (9)-(17)  that 


z** 

-1 

> 

-zt<j) 

-1 

< 

7W 

L\ 

— 

- 

— 

— 

In  consequence  of  (136),  (137),  and  (6),  the  only  non-zero  column 

vectors  IC  and  1^  are  given  by 
n  n  a  j 


ft 

ft 

x-i 

Ji 

i* 

-tt 

-l 

l 

—  — 

(136) 


(137) 


(138) 
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iH 


where  the  column  vector  on  the  right-hand  side  of  (138)  satisfies  (6) 
for  n=l . 

In  view  of  (2)  and  (3),  substitution  of  (138)  into  (1)  and  sub¬ 
sequent  division  by  k  give 


J  r  T.(t)  P.(t) 

—  =  2utcos(D  (I  I  -jg-)  +  2j  u  sinO>  ([  ij  -j—-) 
i  H  |  3  J  j 


(139) 


The  |H  |  written  instead  of  k  on  the  left-hand  side  of  (139)  is  the 

magnitude  of  the  incident  magnetic  field  associated  with  (108).  This 

1 H1 |  is  indeed  equal  to  k.  At  t  =  t  , , ,  the  t  component  of  (139)  re- 
—  P+l 

duces  to 


J  21, 

_t _ lE_ 


|H  I  kp(tp+1) 


cos  (j)  ,  p=l,2, . .  .P-2 


(140) 


At  t  =  t  ,  the  <p  component  of  (139)  reduces  to 
P 

J.  2jlJ 

^-= - sin  <p  ,  p=l,2, . .  .P-1 


In1!  kpp 


(141) 


Here,  Jt  and  are,  respectively,  the  t  and  <J>  components  of  J^.  In  the 

!Jt! 

figures  to  follow,  — —  in  the  <p  -  0°  plane  is  plotted  with  squares  and 

IV  11 ' 

— f—  in  the  <}>  =  90°  plane  is  plotted  with  octagons. 


Figure  4  shows  the  t  and  tp  components 


'V  'V 


and 


of  the  electric 


current  induced  by  the  axially  incident  electric  field  (108)  with  0^=0 


on  an  infinitely  thin  circular  disk  of  radius  0.25X  where  X  is  the  wave- 

lJt!  UJ 

length.  In  Fig.  4,  — — -  is  plotted  with  squares  and  — with  octagons. 

|h‘|  |h‘| 
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Both  quantities  are  plotted  versus  t/A  where  t  is  the  arc  length 
along  the  generating  curve.  The  horizontal  axis  in  Fig.  4  was  labeled 
T/A  because  the  lower  case  letter  t  could  not  be  drawn  by  the  plotter. 

In  Fig.  4,  the  center  of  the  disk  is  at  t  =  0  and  the  edge  at  t  =  0.25A. 
The  electric  currents  in  Fig.  4  and  in  Figs.  5-10  to  follow  were  cal¬ 
culated  with  nfc  =  nT  =  2,  n^  =  20  and  with  the  points  t  ^  ,  j=l,2,...P 
equally  spaced  along  the  generating  curve.  Since  12  octagons  are  in 
Fig.  4,  P=13  therein.  The  electric  current  in  Fig.  4  should  be  twice 
as  large  as  the  magnetic  current  in  Fig.  4  on  page  32  of  [7]. 

Figure  5  shows  the  electric  current  induced  on  a  circular  disk 
of  radius  1.5A  by  the  same  axially  incident  plane  wave  as  in  Fig.  4. 

The  electric  current  in  Fig.  5  should  be  twice  as  large  as  the  mag¬ 
netic  current  in  Fig.  6  on  page  33  of  [7],  Figure  6  shows  the  elec¬ 
tric  current  for  axial  incidence  on  an  infinitely  thin  washer  of 
inner  radius  0.4A  and  outer  radius  1.2A.  The  inner  edge  of  the  washer 
is  at  t  •=  0  and  the  outer  edge  at  t  =  0.8A.  Figure  6  should  be  com¬ 
pared  with  Fig.  3  of  [8].  The  size  of  the  washer  in  Fig.  3  of  [8]  is 
incorrectly  stated.  That  figure  is  actually  a  plot  of  the  electric 
current  on  the  same  washer  as  in  Fig.  6. 

Figures  7  and  8  are  plots  of  the  electric  current  for  axial 
incidence  on  a  cone-sphere  of  cone  angle  20°  and  sphere  radius  0.2A. 
Figure  7  is  for  incidence  on  the  sphere  end  and  Fig.  8  is  for  inci¬ 
dence  on  the  tip  of  the  cone.  The  tip  of  the  cone  is  at  t  *  0.  At 
the  sphere  end,  t  is  approximately  1.48 A.  For  comparison,  see 
Fig.  4.15  on  page  218  of  [9]. 
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Figure  9  shows  the  electric  current  for  axial  incidence  on  an 
open-ended  cylinder  of  radius  A/2tt  and  length  A.  The  plane  wave  is 
incident  on  the  end  of  the  cylinder  for  which  t  *  0.  The  excellent 
results  plotted  in  Fig.  9  here  and  in  Fig.  2.13  on  page  52  of  [2] 
were  both  obtained  by  using  the  electric  field  integral  equation, 
notwithstanding  the  stability  problem  reported  in  [10]. 

Figure  10  is  a  plot  of  the  electric  current  for  axial  inci¬ 
dence  On  the  infinitely  thin  conducting  shell  for  which 

r  *  0.2A 
45°  <  0  <  180° 

where  r  and  0,  being  spherical  coordinates,  are  the  radius  and 
celatitude,  respectively.  This  shell  is  a  spherical  shell  with  an 
axially  symmetric  aperture.  The  pole  of  the  shell  is  at  t  ■  0.  At 
the  edge  of  the  shell,  t  is  approximately  0.471X.  The  plane  wave  is 
Incident  on  the  aperture. 

Numerical  results  for  the  electric  current  on  a  circular  disk 
of  radius  0.02X  not  shown  here  exhibited  a  noticeable  change  in  slope 
near  the  center  of  the  disk.  The  curves  labeled  "a"  in  Figs.  7  and  8 
on  page  34  of  [7]  also  indicate  a  change  in  the  slope  of  the  magnetic 
current  near  the  center  of  the  complementary  aperture.  However,  these 
changes  in  slope  did  not  agree  with  each  other.  Now,  equation  (23) 
of  [11]  does  not  predict  any  noticeable  change  in  the  slope  of  the 
electric  current  near  the  center  of  the  disk  of  radius  0.02X.  The 
changes  in  slope  obtained  by  using  the  computer  program  of  the  present 


liH/r  i 


report  and  the  program  of  [7]  are  obviously  wrong.  The  changes  in 
slope  obtained  by  using  these  programs  are  much  more  pronounced  for 
the  disk  of  radius  0.002A.  However,  they  disappear  when  all  calcu¬ 
lations  are  done  in  double  precision.  Hence,  these  changes  in  slope 

are  due  to  severe  roundoff  error.  This  roundoff  error  occurs  because 

2 

the  vector  potential  terms,  those  containing  the  factor  k  explicit 
in  (90-(12),  are  overshadowed  by  the  rest  of  the  terms  in  (9)-(12), 
the  scalar  potential  terms.  If  these  vector  potential  terms  were  set 
equal  to  zero,  the  moment  matrix  would  be  singular  because  there  are 
several  linear  combinations  of  the  expansion  functions  which  have  no 
electric  charge  associated  with  them. 


PART  TWO 


COMPUTER  PROGRAM 


I .  INTRODUCTION 


The  computer  program  which  Implements  the  numerical  solution 
expounded  in  Part  One  is  described  and  listed  here  in  Part  Two.  This 
program  consists  of  the  subroutine  ZMAT,  the  function  BLOG,  the  sub¬ 
routines  PLANE,  DECOMP,  and  SOLVE,  and  a  main  program.  The  subroutine 
ZMAT  calculates  the  elements  of  the  moment  matrix  in  (6) .  The  func¬ 
tion  BLOG  is  called  by  ZMAT.  The  subroutine  PLANE  calculates  the 
elements  of  the  excitation  vector  in  (6)  for  plane  wave  incidence. 

The  subroutines  DECOMP  and  SOLVE  solve  the  matrix  equation  (6)  for 

IC  and  1^. 
n  n 

The  main  program  obtains  the  electric  current  induced  on  the 

surface  of  the  body  of  revolution  by  the  axially  incident  plane  wave 

(108)  with  8^  -  0  or  Tt  radians.  The  main  program  calls  the  subroutines 

ZMAT,  PLANE,  DECOMP,  and  SOLVE.  It  is  not  difficult  to  generalize  the 

main  program  to  oblique  incidence  because  the  subroutines  ZMAT,  PLANE, 

DECOMP,  and  SOLVE  are  designed  to  calculate  I*"  and  1^  for  n  -  0,1,2,..., 

n  n 

For  the  0-polarized  incident  plane  wave  (108),  1^  is  even  in  n  and  1^  is 
odd  in  n.  For  the  0  polarization  (109),  1^  is  odd  in  n  and  is  even 
in  n.  In  order  to  obtain  far  field  patterns,  the  main  program  must  be 
supplied  with  additional  logic.  This  additional  logic  is  outlined  as 
follows.  According  to  (1-91),  the  far  field  is  obtained  by  premultiplying 
the  solution  vector  to  (6)  by  plane  wave  measurement  matrices  for 
n  *  0,  +1,  +2,...  .  The  plane  wave  measurement  matrices  for  n  *  0,1,2,... 
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can  be  obtained  by  calling  the  subroutine  PLANE.  The  even-odd  behavior 

j«4>r 

in  n  of  the  coefficient  of  e  in  (1-9J)  is  as  follows. 


Receiver 

Polarization 

Transmitter 

Polarization 

Behavior 
in  n 

0 

e 

even  in  n 

<P 

e 

odd  in  n 

e 

$ 

odd  in  n 

<P 

$ 

even  in  n 

. . 

Here,  the  receiver  polarization  denotes  the  component  of  the  far  field 
being  measured.  The  transmitter  polarization  is  the  polarization  of 
the  incident  plane  wave. 


i 
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II.  THE  SUBROUTINE  ZMAT 


The  subroutine  ZMAT(M1,M2,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z)  calcu¬ 
lates  the  moment  matrices  in  (6)  for  n  =  Ml.Ml+l, . . .M2  where  Ml  0 
and  stores  them  in  Z.  Z  is  the  only  output  argument.  The  rest  of 
the  arguments  of  ZMAT  are  input  arguments.  For  n  ■  Ml,  storage  of 

the  Z  submatrices  in  Z  is  as  follows, 
n 

(Z^)^  in  Z(i+N*(j-1)> 

(Z4*),,  in  Z(i+N*(j-1)  +  NP-2) 
n  ij 

(Z^)^  in  Z(i+N*(j-1)  +  (NP-2)*N) 


Here, 


(Z4’4’),.  in  Z(i+N*(j-1)  +  (NP-2)*N+NP-2) 
n  ij 

N  =  2*NP-3 


(142) 


For  n  >  Ml,  the  Z  submatrices  are  stored  in  Z ( (n-Ml) *N*N+1)  to 
n 

Z((n-Ml+1)*N*N)  in  the  same  manner  as  the  Z  submatrices  were  stored 

n 

in  Z(l)  to  Z(N*N)  for  n  =  Ml.  Table  1  relates  the  third  to  eleventh 
arguments  of  ZMAT  to  variables  in  Part  One  of  the  text.  In  Table  1, 
p(t^)  and  zCt^  are  the  values  of  p  and  z  at  t  *  for  i  *  1,2,...P. 


Table  1.  Third  to  eleventh  arguments  of  ZMAT. 


Argument  Variables  in 


of  ZMAT  Part  One 


Minimum  allocations  are  given  by 

COMPLEX  Z(M3*N*N),  G4A(M3),  G5A(M3),  G6A(M3) , 

G4B(M3),  G5B(M3) ,  G6B(M3),  GA(NPHl) ,  GB(NPHI) 
DIMENSION  RH(NP) ,  ZH(NP) ,  X(NPHI) ,  A(NPHI) , 

XT (NT) ,  AT (NT) ,  RS(NP-l),  ZS(NP-l),  D(NP-l) , 
DR(NP-l) ,  DZ(NP-l) ,  DM(NP-l) ,  C2(NPHI), 

C3 (NPHI) ,  R2(NT) ,  Z2(NT),  C4(M3*NPHI),  C5(M3*NPHI), 
C6(M3*NPHI),  Z7 (NT) ,  R7 (NT) ,  Z8(NT),  R8(NT) 
where  M3  -  M2-M1+1. 
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The  elements  of  the  Z  submatrices  are  calculated  according  to 

n 

(48)-(51)  where  G  and  G  ,  are  given  by  (54),  (55),  and  (62)-(66). 
ma  mb 

However,  (62)-(66)  are  modified  through  the  use  of  methods  1,  2,  or  3 
in  the  cases  specified  by  (105)-(107).  The  values  of  ct  and  c^  sug¬ 
gested  in  (104)  enter  via  CT  and  CP  in  lines  10  and  11. 

DO  loop  10  sets 


RS(q)  *  kpq 
ZS  (q)  =  kz 

q 

kA 

D(q)  *  — 2^- 

for  q  ■  1,2,...  NP-1 . 

DO  loop  11  sets 

C2(K)  =  <J>£  and 


DR(q) 

DZ(q) 

DM(q) 


kA 

— sin  v 

2  q 

kA 

— cos  V 
2  q 

A 


20o 


2 

C3(K)  =  4  sin  (y). 


Inner  DO  loop  29  sets 

(n.)  2  V 

C4(M5)  =  it  Aj,  v  sin  (y)  cos(ni{>K) 

_  (n  ) 

C5(M5)  *  9  cos  <}>K  cos(n4>K) 

„  (n.) 

C6(M5)  -  j  Ak  ?  sin  4)r  sin(n^K) 
where  H5  '  K  +  (n-Ml)*NPHI. 

The  calculation  of  (48 ) — ( 51 )  occurs  inside  three  DO  loops  nested 
in  the  following  manner. 


DO  15  JQ  «  1,  MP 
DO  16  IP  -  1,  MP 
DO  31  M  -  1,  M3 
CALCULATION  OF  (48)-(51) 


45 


31  CONTINUE 


16  CONTINUE 
15  CONTINUE 

Here,  JQ,  IP,  and  M  represent,  respectively,  the  variables  q,  p,  and 

(n-Ml+1)  in  (48)-(51).  The  JN  introduced  in  line  52  is  incremented 

in  line  314  so  that  the  subscript  for  (ZtC)  ,  ,  can  be  written  as 

•  n  p-l,q-l 

p  +  JN  when  n  *  Ml.  The  variable  KQ  defined  in  lines  54  to  56  keeps 
track  of  the  cases  for  which  q  =  1  and  q  =  MP.  Because  of  (24),  these 
cases  require  special  treatment.  According  to  (24),  expressions  (48) 
and  (49)  are  absent  when  j  =  q-1  and  q  =  1.  Likewise,  (48)  and  (49) 
are  absent  when  j  =  q  and  q  =  MP. 

The  variables  defined  in  statements  57  to  70  are  needed 
inside  DO  loop  12  or  DO  loop  16.  DO  loop  12  puts  kp'  and  kz'  of  (68) 
and  (69)  in  R2  and  Z2,  respectively.  To  reduce  execution  time, 
references  to  subscripted  variables  as  well  as  calculations  are  being 
done  outside  DO  loops  whenever  possible.  Unfortunately,  this  usually 
increases  the  number  of  statements  and  complicates  the  logic  because 

jk2i  A 

factors  such  as  - q  sin  v  sin  v  in  (48)  are  computed  by  means  of 

a  p  q 

several  statements  scattered  throughout  the  program.  One  way  to  follow 
the  gradual  building  up  of  constants  from  outer  to  inner  DO  loops  is  to 
tabulate  computer  program  variables  versus  variables  in  Part  One  of  the 
text. 

Lines  78  to  88  put  kdQ  of  (102)  in  D6.  Lines  89  to  94  set  KP 
equal  to  the  case  number  in  (103)  and  (105)-(107). 


46 


Lines  96  to  248  put  approximate  values  of  G  of  (56)  in  GmA  for 

ma 

m  =  4,5,6.  Lines  96  to  248  also  put  approximate  values  of  G  ,  of  (57) 

mb 

in  GmB  for  m  =  4,5,6. 

Lines  96  to  174  are  executed  for  case  2  of  (105)  and  for  case  4 

of  (107).  Method  1  is  used  here.  This  method  is  described  by  (71) —(93) . 

The  pure  Gaussian  quadrature  option  for  G  and  G,  advocated  just  after 

a  b 

(93)  is 


G 

a 


n  ,  .  -jkR 

t  (»t) 

tii  v  kRP 


(143) 


n 


t  (O  (n  ) 


-jkR 


Gb  =  £  V* 

£'=1  * 


t  e 


kR 


(144) 


where  R^  has  the  same  meaning  as  in  (81).  In  terms  of  Z7,  R7 ,  Z8,  and 
R8  calculated  by  DO  loop  40, 


kR  =>  Jzl(Z')  +  R7(£’)*(4  sin2 4)) 

p  >  2 


(145) 


kR 


\jz8(.Z')  +  R8U')*(4  sin2(|)) 


(146) 


DO  loop  33  puts  G&  and  G^  in  GA(K)  and  GB(K).  The  index  K  of 

DO  loop  33  corresponds  to  Z  in  ( 91) — ( 93) .  Line  109  puts  k2r2  of  (86) 

P9 

in  RR.  If 


r  +4a 

pq  2  t  q  2  q 


(147) 


then  line  112  sends  execution  to  statement  34  and  G  and  G,  are  calcu- 

a  b 

lated  according  to  (73)  and  (74).  Otherwise,  DO  loop  35  accumulates 

G  and  G.  of  (143)  and  (144)  in  UA  and  UB.  The  purpose  of  the  second 
a  b 
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term  on  the  right-hand  side  of  (147)  is  to  assure  that  the  distance 

between  the  field  point  and  the  closest  point  on  the  line  <J>  =  <J>  on 

the  qth  source  segment  is  no  less  than  —  c^^  before  DO  loop  35  is 

entered.  This  distance  could  be  as  small  as  r  -  4  A  .  DO  loop  37 

pq  2  q 

accumulates  G  ,  and  G,  ,  of  (79)  and  (80)  in  UA  and  UB.  Lines  130  to 
al  bl 

142  add  G  „  and  G,  .  of  (87)  and  (90)  to  UA  and  UB.  Nested  DO  loops 
a  l  bZ 

45  and  46  put  G  of  (91) — (93)  in  GmA  for  m  =  4,5,6.  These  DO  loops 
ma 

also  put  G  ,  in  GmB  for  m  =  4,5,6.  The  DO  loop  indices  M  and  K  cor- 
mb 

respond,  respectively,  to  (n-Ml+1)  and  £  in  (91). 

Lines  176  to  197  apply  method  3  to  G^a>  Expression  (97)  for  G,_a 
consists  of  three  terms,  namely,  two  double  sums  and  a  double  inte¬ 
gral.  Since  the  first  term  in  (97)  is  the  result  of  pure  Gaussian 
quadrature,  the  second  and  third  terms  in  (97)  are  attributed  to 
method  3.  At  this  point,  however,  we  do  not  have  the  first  term  in 
(97),  but  the  modification  of  it  due  to  application  of  method  1  in 
lines  96  to  174.  For  consistency,  the  inner  sum  in  the  second  term  in 
(97)  should  be  replaced  by  the  corresponding  exact  integral  whenever 
(147)  is  true.  This  corresponding  exact  integral  is  given  by 


dt' 


(148) 


Formula  200.01.  of  Dwight  [6]  was  used  to  obtain  the  right-hand  side  of 
(148) .  The  index  K  of  DO  loop  63  corresponds  to  £  in  (97) .  Inner  DO 
loop  65  accumulates  in  D7  the  inner  sum  in  the  second  term  in  (97). 
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Lines  188  and  189  put  (148)  in  D7.  Line  194  puts  in  D8  the  contribu¬ 
tion  due  to  the  second  and  third  terms  in  (97).  DO  loop  67  adds  this 
contribution  to  t he  modified  first  term  in  (97). 

Lines  199  to  248  calculate  G  and  G  ,  according  to  (62),  (63), 

ma  mb 

(64),  (66)  and  either  (65)  or  (94).  The  index  L  of  outer  DO  loop  13 

'jkRnr)t 

corresponds  to  £'  in  (62)  and  (63).  DO  loop  17  puts  (e  v  )/(kR  . 

pit 

in  GA(K)  for  £  =*  K.  If,  in  accordance  with  (106), 

Vq  >  \/(P'  ~  PP>2  +  (Z>  "  ZP>2  (1 

then  (94)  is  used.  Otherwise,  (65)  is  used.  If  (149)  is  not  true, 
then  line  220  sends  execution  to  statement  51.  Otherwise,  lines  221 
to  225  put  in  D6  the  contribution  due  to  the  second  and  third  terms  in 
(94).  Note  that  the  first  term  in  (94)  is  the  right-hand  side  of  (65). 
DO  loop  32  accumulates  (64),  (65),  and  (66)  in  U5,  U6,  and  U7 ,  respec¬ 
tively. 

Inside  DO  loop  31,  lines  262  and  263  put  G^a  and  G^  of  (54)  and 
(55)  in  H4A  and  H4B,  respectively.  Lines  268  to  274  calculate  terms 
in  (48)-(50) .  U5,  U6,  and  U7  belong  in  (48),  U8  and  U9  in  (49),  and 
UC  and  UD  in  (50).  The  variables  K1  to  K8  defined  in  lines  275  to  282 
give  the  locations  in  Z  of  the  matrix  elements  referenced  in  (48)-(50) 
See  Table  2.  In  Table  2,  p  and  q  run  from  1  to  MP  except  where  other¬ 
wise  indicated.  The  forbidden  values  of  p  and  q  in  Table  2  are  due  to 
(23)  and  (24).  The  contributions  (48) —(51)  are  accounted  for  in  lines 
283  to  310.  In  these  lines,  Z(K4),  Z(K6)  and  Z(K8)  are  referenced  for 


Table  2.  Storage  of  matrix  elements  in  Z. 


Location  in  Z 


Z 

Z(K2 
Z(K3 
Z(K4 
Z(K5 
Z(K6 
Z 

Z 


Matrix  Element 


Zn  Vi.q-1 

P 

i  l. 

P 

i  MP 

CVi.s 

P 

i  1, 

ZnC>P.q 

P 

i  MP 

Z**) 

n  p,q-l 

q 

i  1 

z*e) 

n  p,q 

q 

i  MP 

zt<j))  . 

n  p-l»q 

p 

*  1 

Z^) 
n  p,q 

p 

+  MP 

the  first  time,  but  the  rest  of  the  Z's  are  incremented.  The  branch 
statements  interspersed  from  lines  283  to  306  are  due  to  the  for¬ 
bidden  values  of  p  and  q  in  Table  2.  The  seemingly  muddled  and 
repetitive  nature  of  the  Z's  in  lines  283  to  309  is  the  result  of 
an  effort  to  minimize  the  number  of  branch  statements  executed. 


00  lC 
002C 
003 
004 
005 
006 
007 
000 
009 
010 
Otl 
012 
013 
014 
015 
016 
017 
010 
019 
020 
021 
022 
023 
024 
025 
026 
027 
020 
029 
030 
031 
032 
033 
034 
035 
036 
037 
030 
039 
040 
041 
042 
043 
044 
045 
046 
047 
040 
049 
050 
051 
052 
053 
054 
055 
050 
057 
050 
•59 
000 


LISTING  OF  THE  SUBROUT  INC  ZNAT 

THE  SUBROUTINE  ZMAT  CALLS  THE  FUNCTION  BLOC 

SUBROUT INE  ZMATIMI.N2.NP.NPH1.NT.RH.ZH.X.A.XT.AT.Z) 

COMPLEX  Z(l600I.Ul.U2.U3.U4.US.O6.U7.O0.U9.UA.UB.G4At 10).GSA( 1 Ot 
COMPLEX  CMPLX.G6AI I  0) .G4BI 10) >6501  10 1 < GOBI 1 0 1 »M4 A. MS A. Mb A.H4B.HSO 
COMPLEX  H6B.UC. UD. GAl 40 ) • GBI 40) 

DIMENSION  RHI43) .ZHI 43)  .  XI  40  )  .  A<  40  )  .  XT  1 1  0)  •  AT  1 10  )  .RS  <  42  )  .Z5I 421 
DIMENSION  0142 ) • OR (42).DZ( 42 ) • ON  < 4 2) • C 21 40).C3(40) .R21 10)*Z2I 10) 
DIMENSION  C4<200t.C5<200).C6<200).Z7<I6) .R7II0). ZOt I 0)«R0( 10) 
CT*2. 

CP*.  1 

00  10  l*2*NP 
12*1-1 

RSI  I  2  >* . 54| RH| I ) 4RH( 12) ) 

ZSI 12)*. S«t ZHI I ) 4 ZHI 12)1 
Ol*.S4<  RHI I )-RH| 12)) 

02*. 5*1 ZHI I )— ZHI 12)1 
01  I  2 )=SCRTI0|40 1 >02402) 

ORI I2)*0t 
0Z1 1 2 )*02 

0M| 12 )*0 I 121/RSI 12) 

10  CONTINUE 
M3*M2-MI4| 

M4*MI— 1 
PI  2* 1.570790 
OO  It  K* 1 . NPHl 
PHsP I 24| X1K)4|.) 

C2CK )*PH4PH 
SN*SINI.S4PH) 

C3(K)*4.4SN*SN 
AI*PI24AIK) 

04*.S4A1 4C3CK) 

05-AI4C0SIPH) 

06*A 1 4SI NIPH) 

M5*K 

OO  29  M*1.M3 
PHM*(M44M)4PH 
A2*C0S I PHN) 

C4 1  MS 1*044 A2 
CS1M5 )=DS4A2 
C6I M5)*064SINI PHM) 

MS*MS4NPHI 
29  CONTINUE 

11  CONTINUE 
MPsNP-1 
m  r*MP-i 
N*NT4NP 
N2N*NT4N 
N2*N4N 
U1*(0...S) 

U 2*10. .2.) 

JN»- I— N 
OO  15  J0*1.MP 
K0*2 

IFIJO.EQ.I)  XO*l 
IF I JO.EO.MP)  KQ*3 
R|*RSI JO) 

Z1*ZS< JO) 

01*01  JO) 

02*DRI JO) 
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061 

062 

063 

064 

06S 

066 

067 

060 

069 

070 

071 

072 

073 

074 

07S 

076 

077 

070 

079 

000 

06t 

002 

003 

004 

OOS 

006 

007 

oea 

009 

090 

091 

092 

093 

094 

095 

096 

097 

090 

099 

100 

101 

102 

103 

104 

105 

106 
107 

too 

109 

110 
111 
112 

113 

114 

115 

116 
117 
110 

119 

120 


1 

i 


03*02(301 
04*02/01 
05*001301 
SV*02/0t 
CV*OS/01 
T6*CT*0t 
T62*T6»01 
762*762*702 
R6*CP*RI 
062*06*06 
00  12  L«1«NT 
02(1.1*01  *02*  271  LI 
Z2  <L)-Z1  403*XT(0 
12  CON7INUE 
U 3-02*01 
U4-D3*UI 
00  16  (P*UNP 
R3*R$(IP> 

Z3*ZS(1P> 

R4— R1-*R3 
Z4*2l-23 
FM-R**SV4Z4*CV 
PHM— A8S( FM1 
PH-46SI R4*CV-Z4*SV> 

06— PH 

IF(PHM.LE.Dl)  60  TO  20 
D6*PHM>0I 

06—S0R71 06*064PH*PM) 

26  (FKP.eO.JOl  60  TO  27 
KP*t 

IF176.GT.06)  HP* 2 
IF1R6.6T.06I  KP*3 
GO  70  20 

27  KP*4 

20  GO  TO  ( 4  I *42*41 «42) •OP 
42  DO  40  L*  !•  NT 
D7-R21L1-R3 
D8-Z21LI-Z3 
Z7 (L 1 *07*07400*00 
R7(L»*R3*R2U.I 
Z8(L  1* •25*2711. ) 

R8CL1-.  25*070.1 
40  CONTINUE 

Z4*R4*R4 424*24 
R4*R3*R1 
R5«.S*R3*SV 
00  33  K»  l.NPHI 
A1*C3(K) 

RR-244R49AI 
UA*0  • 

uo*o* 

IFtRR.LT.T62)  60  TO  3* 

DO  33  L»1 .NT 

R*S0R7( Z  7(L )4R7 (L)*At ) 

SN—S1N1R1 

CS*COS(R) 

UC*AT(L )/R*CMPCX(CS« SNl 
UA-UA4UC 
UB*XT(L)*UC4UB 
35  CONTINUE 
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121 

00  TO  36 

122 

34 

00  37  UltNT 

123 

R=SORT ( Z8(L )*R8(L)4A I ) 

124 

SN*— SINCRI 

125 

CS*COS(R) 

125 

UC*AT(L>/R4SN4CMPLX(-SN.CSI 

127 

UAsIMMJC 

128 

UB*XT<L)4UC*U8 

129 

37 

CONTINUE 

130 

A2*FM*R54A1 

131 

09*RR-A2*A2 

132 

R-ABS(A2> 

133 

D7*R-0l 

134 

DB*R»Ot 

135 

D6*S0RT< 08400*09) 

1 34 

RsSORT ( 07407*09  ) 

137 

IF(D7.GE.O.)  GO  TO  38 

138 

Al*At.OG(  (08*06  )  *(—07*R  1/091/01 

139 

GO  TO  39 

140 

38 

A|*A10G( (00*06) /(07*Rt)/0t 

141 

39 

UA*At*UA 

142 

U  B*  A  2  *  (  4  ./ 1  06*R  )  -  A 1 )  /D 1  *  UO 

143 

36 

GA(K)*UA 

144 

G8(K)*UB 

145 

33 

CONTINUE 

146 

Kl*0 

147 

OO  45  N*I.M3 

148 

H4A*0. 

149 

H5A*0. 

ISO 

H6A*0. 

151 

H4B»0. 

152 

H58-0* 

153 

H6B*0. 

154 

OO  46  K-l.NPMI 

155 

81*8141 

156 

06*C4CK1) 

157 

07-C5181 > 

158 

08*C6(Kt I 

159 

UA*GA(K) 

160 

U8*G8(K» 

161 

M4A*064UA*H4A 

162 

H5A*074UA*HSA 

163 

H6A«0e*UA*H6A 

164 

M4 8*  064  US*M4B 

165 

HS83074UB*HS8 

166 

H6B*084UB*H66 

167 

46 

CONTINUE 

168 

G4A(M)*H4A 

169 

G5A(M)*H5A 

170 

G6A(M)*H6A 

171 

G4B(M)*H48 

172 

GSB(M)*HS8 

173 

G6B( M)*H6B 

174 

45 

CONTINUE 

ITS 

IFfXR.NE.4)  GO  TO  47 

176 

A2*0t/(Pt24flil 

177 

06*2. /Ol 

178 

08*0. 

179 

OO  63  K* 1 tNPMl 

180 

A1*R44C2(R» 

*  J  n'SilfV  i  liliml 
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(•1 

R*R4*C3(K) 

182 

1FIR.LT.T62)  60  TO  66 

182 

07*0. 

18* 

OO  65  L* l.NT 

189 

D7*D7*AT (L l/SQRTl 2?(L)*A 1 ) 

186 

69 

CONTINUE 

187 

GO  TO  66 

188 

66 

AI*A2/(  X(K)*1. t 

189 

07*06 *ALOG( Al *SQRT( 1 .*AI6AIII 

190 

66 

0S*0S*A(K»*D7 

191 

63 

CONTINUE 

192 

AI*.5*A2 

193 

A 2*1 ./At 

196 

08*— P  f 2*09*2./RI*(6L0Gt A2 )*A2*BL0G(AI ) ) 

(99 

00  67  M* I .M3 

196 

G9A(M)*D8*GSA(N) 

197 

67 

CONTINUE 

198 

GO  TO  67 

199 

61 

OO  29  M*1.M3 

200 

G6AI M)*0. 

201 

GOA  (10*0* 

202 

G6A( M)*0. 

203 

G46(M)*0. 

206 

GS8(M)*0. 

209 

G69(M)*0. 

206 

25 

CONTINUE 

207 

OO  13  L*l.NT 

208 

A1*R2(L) 

289 

R6*A 1-R3 

210 

Z  6*22(0-23 

211 

Z4*R4*R4*246Z4 

212 

R4*R36A1 

2(3 

OO  17  K* 1 .NPHI 

216 

R*SORT( Z4*R4*C3(K) ) 

219 

SN*-S(N(R) 

216 

CS*COS( R) 

217 

GA(K )*CMPLX(CS« SNl/R 

218 

17 

CONTINUE 

219 

06*0. 

220 

IF (R62.LE.Z6)  GO  TO  St 

221 

OO  62  K*1.NPHI 

222 

06*06* A ( K I/SORT ( Z4*R4*C2(K ) ) 

223 

62 

CONT I NUE 

226 

Z4*3 . 141 S93/SOR  T( Z6/R4) 

225 

06*— R  12*06*  ALOGI  Z4*$QRT(  l  •  *24*24  DFS0RTIR4) 

226 

91 

A  t*AT(L) 

227 

A2«XT(LI*A1 

228 

Kt*0 

229 

OO  30  M*1«H3 

230 

U9*0. 

231 

U6*0. 

232 

U7*0. 

233 

OO  32  X*t .NPHI 

236 

UA*GA(K) 

239 

K1*K1*1 

236 

U5*C4(K1 I6UA6US 

237 

U6*C5(K1 >*UA*U6 

238 

U7*C6(KI I6UA6U7 

239 

32 

CONTINUE 

260 

U6*06*U6 
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241 

G4A( M>*AI*US7G4A(M| 

242 

654<  N 1*4 1 *U67G5A< Ml 

243 

G6A(M>*At*UT7G6A(MI 

244 

GABf  M|*A24US 7048(8} 

245 

G5atM>*A2*U67G5BtM) 

246 

G68( M>*A2*U77G68(M» 

247 

30  CONTINUE 

246 

13  CONTINUE 

249 

47  AI*OP(IP» 

250 

UA*AI«U3 

251 

UB*OZ(f PI9U4 

252 

A2*Of IPI 

253 

06*- 42*02 

254 

07*0 I4AI 

255 

08*01 *42 

256 

257 

00  31  M*I«MJ 

256 

FM*M47N 

259 

AI*FM«0M( IP) 

260 

HS4*GSA(  Ml 

261 

HSB*GS8( Ml 

262 

H4 A*C4A( MI4H56 

263 

H4B*G4B(M|7HS8 

264 

H6A*G6A( Ml 

265 

H68*G68f  Ml 

266 

U7*UA*H5A4U84N4A 

267 

U8*UA*HS87U86H4S 

266 

U5*U7-U8 

269 

U6*U77U8 

270 

U7*-UI*H4A 

271 

U8*064H64 

272 

U9-064H68-AI4H4A 

273 

UC*07«( H6A7044H6S) 

274 

UO*FM4DS4H4A 

275 

Kt*(P7JM 

276 

K2*R|4| 

277 

K3*K 1 7N 

278 

K4*K27N 

279 

KS*K27NT 

280 

K6*K4«NT 

281 

K7*K37N2N 

282 

K6*K47N2N 

283 

GO  TO  (I8.20.I9I.KO 

284 

18  Z(K6)*U8»U9 

301 

285 

IFIIP.EQ.II  GO  TO  21 

302 

286 

Zl K3)*Z( K3I 7U6-U7 

303 

287 

Z( K7I*Z ( K7| 4UC-U0 

304 

288 

IFf (P.EQ.MPI  GO  TO  22 

305 

289 

21  Z(K4)*U67U7 

306 

290 

Z(K8l*UC7UO 

307 

291 

60  TO  22 

308 

292 

19  Zf K5I*Z( KSl 706-09 

309 

293 

tFflP.EO.il  GO  TO  23 

310 

294 

ZfKI }*Z(K1 I 7U57U7 

311 

295 

Z(K7l*Zf K7I7UC-U0 

111 

296 

IFf (P.EQ.MPI  GO  TO  22 

313 

297 

23  ZIK2I*Z(K2|7US-U7 

314 

296 

Z(KSI*UC7UO  ’ 

319 

299 

GO  TO  22 

316 

390 

20  Z(KSI*Z(K5I«U6-U9 

317 

Z(Kt)*U8«U9 
(FflP.EQ.il  CO  TO  M 
ZCK»)aZ(M»*W*ur 
Z(K3l*Zt  K3I7U6-U7 
Z t K7 }*Z ( KT I 7UC-U0 
IFf (P.EQ.MPI  CO  TO  IS 
Z*  Z(«>»Z(«J*U5-W 
Z(K4I*U67U7 
Zf K8)*UC7U0 

22  ZIM»MT).U2»CM»(HSMO«H58l-*l*«>l 
JM*JN7N2 
31  CONTINUE 
16  CONTINUE 
JN*.JN7N 
IS  CONTINUE 
RETURN 
END 
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III.  THE  FUNCTION  BLOG 

The  function  BLOG(x)  calculates  log  (x  +  /I  +  x2)  for  x  >_  0. 

If  x  is  appreciable  compared  to  1,  the  FORTRAN  supplied  subroutine  for 
the  logarithm  suffices.  However,  if  x  is  much  smaller  than  1,  this 
subroutine  fails  because  of  excessive  roundoff  error.  From  formulas 
700.1.  and  706.  of  Dwight  [6], 

log(x  +  /£+**)  «  x(l  -  x2  +  2^75  x*  -  2-i»i?7  **  +  *  *  •> »  *2  <  1  <150> 

If  |x|  £  .1,  the  approximation 

_ _  2  4 

log (x  +  /l+x2)  -  x(l  -  (151) 

incurs  an  error  of  less  than  one  part  in  10^.  The  function  BLOG(x) 
uses  the  FORTRAN  supplied  subroutine  for  the  logarithm  for  x  >  .1  and 
(151)  for  x  <  .1. 


1=  LISTING  OF  TtC  FUNCTION  BLOG 

2  FUNCTION  BLOGIXt 

3  IFCX.GT..1I  GO  TO  I 

*  X2«X*X 

5  BLOG* Cl  .07S*X2-,  1666667)  *X2M.)*X 

6  RETURN 

7  t  BLCG*ALOG(X*SQRT <I.+X*X>) 

RETURN 
ENO 
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IV.  THE  SUBROUTINE  PLANE 


The  subroutine  PLANE(M1,M2,NF,NP,NT,RH,ZH,XT,AT,THR,R)  calculates 
the  elements  of  the  plane  wave  excitation  vectors  according  to  (124)-(127) 
and  (132)-(133)  and  stores  them  in  R.  R  is  the  only  output  argument. 

The  rest  of  the  arguments  of  PLANE  are  input  arguments.  There  are  NF 
angles  of  incidence  0  of  (110)  and  n  -  Ml,  M1+1,...M2  where  Ml  _>  0. 

The  Kth  angle  of  incidence  resides  in  THR(K)  in  radians.  For  the  first 
angle  of  incidence  and  for  n  *  Ml,  storage  in  R  is  as  follows. 


Here, 


in  R(i) 

-V^®  in  R(i+NP-2) 
-V^  in  R(i+N) 

Vnt  in 

N  =  2*NP-3 


(152) 


(153) 


The  minus  signs  are  attached  to  and  in  (152)  so  that,  accord- 

ni  ni 

ing  to  (1-100)  and  (1-104),  the  vectors  stored  in  R  will  be  measurement 
vectors.  For  the  Kth  angle  of  incidence  and  for  n  _>  Ml,  the  storage 
arrangement  of  V*®,  -V^ ,  -V^,  and  is  still  the  same  as  indicated 
above,  but  the  storage  area  now  extends  from  R(2*N*((K-1)*(M2-M1+1)  + 


n-Ml)  +  1)  to  R(2*N*((K-1)*(M2-M1+1)  +  n-Ml+1))  instead  of  from  R(l)  to 
R(2*N).  Table  3  relates  the  fourth  to  ninth  arguments  of  PLANE  to 
variables  in  Part  One  of  the  text.  In  Table  3,  p(t^)  and  z(t^)  are  the 
values  of  p  and  z  at  t  *  t^  for  i  -  1,2,...P. 
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Table  3.  Fourth  to  ninth  arguments  of  PLANE 


Argument  Variable  in 

of  PLANE  Part  One 


Minimum  allocations  are  given  by 

COMPLEX  R(2*N*NF*(M2-M1+1)),  FA(M2+3) ,  FB(M2+3) 

DIMENSION  RH(NP) ,  ZH(NP),  XT (NT) ,  AT (NT) , 

THR(NF),  CS(NF) ,  SN(NF) ,  R2(NT),  Z2(NT) 

where  N  is  given  by  (153). 

The  index  IP  of  DO  loop  12  obtains  p  in  (124)-(127).  DO  loop  13 
puts  y  of  (134)  and  kz^  of  (135)  in  R2(L)  and  Z2(L),  respectively, 
for  l  ■  L.  The  index  K  of  DO  loop  14  obtains  the  Kth  angle  of  inci¬ 
dence. 

The  index  L  of  DO  loop  15  obtains  l  in  (132)  and  (133).  Line  48 
puts  2  sin  in  X.  Lines  49  to  73  calculate  S  and  BJ(m+2)  so  that 


58 


BJ(m+2)  *  S*J  (k3„  sin  0  ),  m  -  Ml-1,  Ml,  ...  M2+1 
m  %  t 

m  ?  -1 

If  the  argument  of  the  Bessel  function  J  in  the  above  equation  does 

m 

not  exceed  10  ^ ,  lines  50  to  54  use  the  approximations 


m  ^  0 


1  ,  m  =  0 


in  order  to  obtain  BJ(m+2)  and  S.  The  purpose  of  lines  56  and  57  is  to 

—8 

obtain  M  so  large  that  |JU  «(k0„  sin  0)1  is  roughly  10  .  Line  58 

x.  t 

assures  that  M  is  at  least  as  large  as  M2+3.  Lines  59  to  67  start  with 


JM-2(x)  =  0 

jm_3(x>  =  1 


and  use  the  recurrence  relation 


Vl<x>  '  T  V*’  '  Jn+l<x) 


taken  from  (9.1.27)  on  page  361  of  [12]  to  calculate  Jr(x)  for  n  *  M-4, 
M-5,...  0.  Lines  68  to  73  use 

1  -  JQ(x)  +  2J2(x)  +  2J4(x)  +  2J6(x)  +  ... 

taken  from  (9.1.46)  on  page  361  of  [12]  to  obtain  the  normalization 


constant  S.  As  the  index  of  DO  loop  15  changes,  DO  loop  25  accumulates 

F  and  F  .  of  (132)  and  (133)  In  FA(m+2)  and  FB(m+2),  respectively.  If 
ma  mb 

F  ,  and  F  ,  .  are  needed,  lines  83  and  84  use  the  formulas 

“1  yd  f  D 
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in  FA(1)  and  FB(1),  respectively. 


to  store  F  ,  and  F  ,  , 

-l,a  -1 , b 

With  reference  to  (124)—(127) ,  the  index  M  of  DO  loop  27  obtains 
(n+2).  Inside  DO  loop  27,  UA  is  tt j n .  The  variables  U2,  U3,  U4,  and  U5 


calculated  in  lines  95  to  98  are  needed  in  order  to  assemble  the  right- 


hand  sides  of  (124)  and  (126).  The  variables  Kl,  K2,  K4,  and  K5  are  the 

subscripts  of  R  for  .  and  respectively.  Lines 

n,p-l  np  n,p-l  np 

102  and  103  obtain  (125)  and  (127).  The  branch  statement  in  line  104  is 

^  0  t  (t) 

necessary  because  neither  V  ,  nor  Vc?  ,  exists  for  p=l.  In  lines 

n,p-l  n,p-l 

105  and  106,  ,  and  Vt^>  ,  are  incremented.  The  branch  statement  in 

n, p-1  n,p-l 

line  107  is  necessary  because  neither  nor  Vt(^  exists  for  p  =  NP-1 .  In 

np  np 

lines  108  and  109,  and  Vt<®>  are  referenced  for  the  first  time. 

np  np 


001 c  LISTING  OF  THE  SUBROUTINE  PLANE 

002  SUBROUTINE  PLANE! Ml . M2.NF. NP.NT. RH.2H. XT. AT.THA.lt) 

003  COMPLEX  RI240I .U.Ul .UA.UB.FA! 1 0} .FBI 1 0  I .F2A.F28.FI A.  FIB. U2.U3.U4 

004  COMPLEX  U5.CMPLX 

005  01  PENSION  RHI43I .ZH! 43 ) . XT ! 1 0 > . ATI 10 1 . THft 1 3) .CS! 3) «SN| 3)  .ft  2!  101 

00b  DIMENSION  221 1  0) .BJ( SOI 

00 7  MP*NP-I 

008  NT-MP-I 

004  N*MT4MP 

010  N2*24N 

Oil  DO  I  1  K* l.NF 

012  X=THR|K| 

013  CS ( K ) =C0S(  XI 

014  SNIK)*SINfX) 

OtS  tl  CONTINUE 
016  U*I0  .  .1  « I 

017  UI*3.141S93«U**N1 

018  N3>MI»| 

019  N4*M243 

020  IFCMt.EO.OI  M3  a  2 

021  NS*MI »2 

022  M6=M2*2 

023  DO  12  IP*1.MP 

024  K2*(P 

025  I=IP*1 

026  OR*. 5*<RH( I >-RHl IP)) 

027  OZ*.S*IZH|I)-ZH|IP)» 

028  Ol*S0RT(0R40R4DZ*OZt 

029  R1*.2S*1RH1  IH-RH11P)) 

030  Zl=.S41ZHlll4ZHltPII 

031  DR*. SADR 

032  D2=DR/RI 

033  DO  13  L*1.NT 

034  R21L1=RI ♦DR4XTILI 

035  22<L»*Z|40Z4XT(L> 

036  13  CONTINUE 

037  DO  14  K*t.NF 

038  CC-CSOCI 

039  SS*SNIKI 

040  03*DR4CC 

041  04*— 0Z4SS 

042  D5*01 4CC 

043  DO  23  M*N3.N4 

044  FA!NI*0. 

045  FB! M1*0* 

046  23  CONTINUE 

047  DO  15  L* I .NT 

048  X*SS«R2(L) 

049  IFIX.GT..5E-7>  GO  TO  19 

050  DO  2  0  M*M3.M4 

051  BJ1 M 1*0. 

052  20  CONTINUE 

053  83! 2 1*1 . 

054  5*1. 

055  GO  TO  IS 

056  19  M*2.84X4l4.-2./X 

057  IF ! X.LT. .5)  M*l 1 .84AL0GI  0! X) 

058  IF1M.LT.M4)  M-M4 

059  8  JIM) *0. 

060  JM*M- I 


061 

BJI JNt«l . 

062 

DO  16  J*4.M 

063 

J2*JM 

064 

JN.JM-l 

06S 

Jt-JM-l 

066 

BJ< JM)»Jl/X*BJC J2I-BJC JMF2I 

067 

16 

continue 

068 

s*o. 

064 

IFtM.LE.4t  GO  TO  24 

070 

OO  17  JS4.M.2 

07| 

S«  S4-BJC.lt 

072 

17 

CONTINUE 

073 

24 

S^BJC  2 14-2.65 

074 

18 

ARG«Z2<L1»CC 

075 

UA^ATCL)  /S4CMPLX  (COS  1  ARG I  *  S  INI  Aft  G)  1 

076 

UB=XT(Lt*UA 

077 

DO  25  MsM3.M4 

078 

F  AC  Mt  *B J ( Mt 6UA4F AIM) 

079 

FBI M l=BJ ( Mt  6UB4FBIM1 

080 

2S 

CONTINUE 

081 

IS 

CONTINUE 

oe2 

tFIMI.NE.Ot  CO  TO  26 

083 

FA(tl=FA<31 

084 

FB( 1 1*FBC3I 

OSS 

26 

UAaUt 

086 

DO  27  M>MS.M6 

087 

M7»l»-t 

088 

M8*M4| 

089 

F2  A»UA*  (FAC  MS  >4FAI M7 1 1 

090 

F2B=UA*(FBt M81AFBIM71 1 

091 

UB»U*UA 

092 

F  t A»UB* <  FAI MSI— FACM7) t 

093 

F 1 8*UB*< FBI M8I-F8CM7) » 

094 

U4«04*UA 

095 

U2SD34F t A4U44FAIM) 

096 

U3*D3*F t  B4U44F8I M> 

097 

U4»DR4F2A 

098 

U5=0R*F2B 

099 

KI*K2-t 

too 

K4*K14N 

tot 

K5SK2  4N 

102 

RC  X24-MT )=-OS*(F2A402*F2Bl 

103 

R(K54MT)>DI*tFIA»02»FIB) 

I  04 

(Ft  IP.EO.lt  CO  TO  21 

tos 

RtKt  >*RCKll4-U2-U3 

106 

R  |  K4 1  »R  (  K4 1 4-U4-U5 

107 

IFIIP.EQ.MPt  GO  TO  22 

too 

2t 

R(K21«U2HI3 

109 

R(K5t«U44US 

tto 

22 

K2*K24-N2 

t  tl 

UA«UB 

t  12 

27 

CO  NT  (NUE 

1  13 

14 

CONTINUE 

114 

12 

CONTINUE 

tts 

RETURN 

116 

END 
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V.  THE  SUBROUTINES  DECOMP  AND  SOLVE 

The  subroutines  DECOMP (N,  IPS,  UL)  and  SOLVE (N,  IPS,  UL,  B,  X) 
solve  a  system  of  N  linear  equations  in  N  unknowns.  The  input  to 
DECOMP  consists  of  N  and  the  N  by  N  matrix  of  coefficients  on  the  left- 
hand  side  of  the  matrix  equation  stored  by  columns  in  UL.  The  output 
from  DECOMP  is  IPS  and  UL.  This  output  is  fed  into  SOLVE.  The  rest  of 
the  input  to  SOLVE  consists  of  N  and  the  column  of  coefficients  on  the 
right-hand  side  of  the  matrix  equation  stored  in  B.  SOLVE  puts  the 
solution  to  the  matrix  equation  in  X. 

Minimum  allocations  are  given  by 
COMPLEX  UL(N*N) 

DIMENSION  SCL(N) ,  IPS(N) 
in  DECOMP  and  by 

COMPLEX  UL(N*N) ,  B(N) ,  X(N) 

DIMENSION  IPS (N) 

in  SOLVE. 

More  detail  concerning  DECOMP  and  SOLVE  is  on  pages  46-49  of  [13]. 


001  C 

LISTING  OF  THE  SUBROUTINES  OECOMR 

AND 

SOLVE 

•  02 

SUBROUTINE  OECOMPIN. IPS* UL 1 

003 

COMPLEX  UL<I600I<PIV0T<EW 

00* 

DIMENSION  SCL<401<1PS<401 

005 

00  5  1*1 <M 

006 

!PS<  11*1 

007 

RN*0* 

ooa 

31*1 

009 

DO  2  J*  1  <N 

010 

ULM*ABSt  RE  ALIULI  Jill  l»ASS  <AIMAG<UL<UD1  > 

Otl 

JI*JI«N 

012 

(FIRN-ULM1  1.2*2 

0t3 

1  RN*ULN 

01* 

2  CONTINUE 

015 

SCL< 11*1 */RN 

016 

5  CONTINUE 

017 

NNl*N-t 

016 

1(2*0 

019 

DO  17  K*1.NMI 

020 

BIG*0* 

021 

OO  11  I*K.N 

022 

1P*IPS(1 1 

023 

IPK*IP*K2 

02* 

S 1 2E* 1 ASS<  REAL I UL I  IPX 1 1 1 » ABS  <  A I MAGI UL< IPK 1 1 1 1 *SCL< IP 1 

025 

IFIS1ZE-BIG)  11.11*10 

026 

10  BIG*S1ZE 

0  27 

1PV-I 

026 

11  CONTINUE 

029 

IFI1PV-K1  1*. 15.1* 

030 

t*  J* IPSIK 1 

031 

|PS(KI*(PS< IPV1 

032 

IPS! (PV1*J 

033 

IS  KPP*1PS<K1*K2 

03* 

PI VOT*UL 1KPP1 

035 

KPt*K*l 

036 

DO  16  t*KPl.N 

037 

KP*KPP 

038 

IP*|PS<I 1*K2 

039 

EM*- ULl IP 1/PIVOT 

0*0 

IB  ULl IPf*—EN 

0*1 

DO  16  J*XP1,N 

061 

DO  1  J*|,!N1 

0*2 

IP*|P*N 

062 

SUM*SUM*UL<IP|*X<J» 

0*3 

KP*KP*N 

063 

1 

I P* I P*N 

0*4 

ULl 1PI*UL< IP l+ENPULl KP1 

06* 

2 

XI I1*BI IPB1—SUM 

0*5 

16  CONTINUE 

065 

K2*N*<N— 1 1 

0*6 

K2*K2*N 

066 

IP*I  PS(  NXK2 

0*7 

17  CONTINUE 

067 

XI N1*X<  N l/ULI IPl 

0*6 

RETURN 

066 

DO  •  f BACK*2*N 

0*9 

END 

069 

t*NP|— IBACK 

050 

SUBROUTINE  SOLVE<N. tPS.UL* B.XI 

070 

K2*K2-N 

051 

COMPLEX  UL 1 16001 *6(401 *X< 40! .SUM 

071 

IPI* IPS< I |*K2 

052 

DIMENSION  IPS1401 

072 

IP  1*1*1 

053 

NP 1*N*1 

073 

SUM*0. 

OS* 

IP*IPS<II 

07* 

IP*IPI 

OSS 

XI ll-BIIPI 

075 

00  3  J* IPt *N 

056 

00  2  1*2. N 

076 

IP«tP»N 

057 

IP*IPS<II 

077 

3 

SUM»SUM*UL < IPl • X 1 J1 

OSS 

IPS* IP 

076 

* 

Xlt)*(X< II-SUMI/ULIIPII 

059 

IMI*I-I 

079 

RETURN 

060 

SUM*0* 

060 

END 
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VI.  THE  MAIN  PROGRAM 


The  main  program  calculates  the  electric  current  induced  by  a 

plane  wave  axially  incident  on  a  perfectly  conducting  surface  of 

revolution.  This  plane  wave  is  given  by  (108)  with  0  *  0  or  it  radians. 

The  components  of  the  electric  current  are  obtained  from  (140)  and  (141) 

in  which  if  and  if  are  the  pth  elements  of  the  vectors  if  and  if  which 
Ip  lp  11 

satisfy  (6)  for  n=l. 

Punched  card  data  are  read  in  according  to 
READ(1,15)  NT,  NPHI 

15  FORMAT (2 13) 

READ(1,10) (XT(K) ,  K-l,  NT) 

READ(1,10) (AT(K) ,  K-l,  NT) 

10  FORMAT (5E14. 7) 

READ (1,10) (X(K) ,  K-l,  NPHI) 

READ(1 , 10) (A(K) ,  K-l,  NPHI) 

READ (1,16)  NP,  BK,  THR(l) 

16  FORMAT (13,  2E14.7) 

READ (1,18) (RH(I) ,  1=1,  NP) 

READ(1,18) (ZH(I) ,  1=1,  NP) 

18  FORMAT (10F8. 4) 

Here,  BK  is  the  propagation  constant  k  ah'd  THR(4-X-~is_...the  angle  of 
incidence  in  radians.  THR(l)  must  be  either  0  or  it.  The  input 
variables  NT,  NPHI,  XT,  AT,  X,  A,  and  NP  are  defined  in  Table  1. 

These  input  variables  can  therefore  be  fed  directly  into  the  subroutine 
ZMAT.  However,  RH  and  ZH  must  be  multiplied  by  BK  before  being  fed 


into  ZMAT.  More  precisely,  RH  and  ZH  are  values  of  p  and  z  so  that 
the  product  of  RH  with  BK  is  the  RH  in  Table  1,  and  the  product  of  ZH 
with  BK  is  the  ZH  in  Table  1.  The  sample  input  and  output  data  listed 
along  with  the  main  program  are  for  the  spherical  shell  of  Fig.  10. 

Minimum  allocations  are  given  by 

COMPLEX  Z(N*N),  R(2*N) ,  B(N),  C(N) 

DIMENSION  RH(NP) ,  ZH(NP),  X(NPHI), 

A(NPHI) ,  XT (NT) ,  AT (NT) ,  IPS(N) 

where  N  =  2*NP-3. 

With  reference  to  (6),  line  41  puts  the  moment  matrix  in  Z. 

Line  46  puts  the  excitation  vector  of  (6)  and  the  negative  of  the 
excitation  vector  of  (6)  in  R(l)  to  R(2*NP-3).  These  excitation 
vectors  are  for  the  0-polarized  plane  wave  (108)  and  their  elements 
are  called  and  Storage  in  R  is  according  to  (152).  Now, 

-V^  and  are  also  stored  in  R,  but  are  not  used.  Lines  47  to  52 
put  and  in  B.  Lines  55  and  56  put  the  solution  vectors  I*  and 
1^  to  (6)  in  C.  DO  loop  24  prints  out  (140)  at  <{>  ■  0®.  DO  loop  27 
prints  out  (141)  at  <J>  *  90°. 


001  C  LISTING  OF  THE  MAIN  PROGRAM 

002  C  THE  SUBROUTINES  2MAT •  PLANE.  OECOMP.  AND  SOLVE  ARE  CALLED* 

003  / /PGM  JOB  I XXXX. XXXX . 1*2).*MAJVZ .JOE*. REGION* ZOOK 
004//  EXEC  WATFIV 
OOS//GO.SVSIN  00  4 

006  *  JOB  MAUTZ.T1ME*3.PAGES*60 

007  COMPLEX  Z( I  6001 *RI 2401 .B I  4  01 >C{ 40) .U.C I 

006  DIMENSION  THRO)  .RHI 43)  •  ZHI 43 1.XI46).  A  (40).  XT  I  101. ATI  I0).|PS(40» 

009  REAOII.IS)  NT.NPHI 

010  IS  FORMAT 12 131 

Oil  WRITEI3.30)  NT.NPHI 

012  30  FORMAT! •  NT  NPHI •/! X. 13. (SI 

013  RE AO II.IOIIXTIX). K*l .NT) 

014  READ! 1.10)1 ATIK ) .X* I .NT) 

OIS  10  FORMAT I 5EI 4*7) 

016  VRITE13.  IDIXTC  K).**I.NT) 

017  WRITE  13  » 12 ) t  AT IK). K* I .NT ) 

010  II  FORMAT! •  XT'/!  1 X.SE14.7) ) 

019  12  FORMAT!*  AT  *  fl I X. S£ 14.7) ) 

020  RE  AO  I  1  •  1  0)  I XI K)  •  X»1  •  NPHI  I 

021  READ! t. 101! A(K) «K» I • NPHI  I 

022  WRITEI3. I3IIXIK) «K»1 .NPHI I 

023  WRITE!3.l4)!A!K).Ks|.NPHC> 

024  13  FORMAT! •  X* /! I X.SEI4 *7*1 

02S  14  FORMAT!'  A'  /I 1 X. SE14.7) ) 

026  READ! 1.16)  NP.BK.THRII) 

027  16  FORMAT! I 3.2E14.7) 

026  WRITEI3.I7)  NP.SK.THRIt) 

029  17  FORMAT!'  NP* .6X.* SK' . 12X. *T HR •/ IX. 13. 2EI4.7) 

030  RE AO I  1. 1811 RHI I ) . 1*1 .NP) 

031  READ! I . 1 B) I ZHI I )■ 1*1 .NP) 

032  IS  FORMAT! I 0F8.4) 

033  WRITEI3. 19IIRH! (>.t*t.NP) 

034  WRITE  13.20 ) IZH! t)«(*l*NP) 

035  19  FORMAT!*  RH*/{ IX. I0F0.4) I 

036  20  FORMAT!*  ZH* /I 1 X • 10F8.4) ) 

037  DO  28  J=I*NP 

036  RHI  J  )*BK4RH|  J) 

039  ZHI J)*BK»ZH«J> 

040  20  CONTINUE 

041  CALL  ZMATI t . I .NP. NPHI .NT .RH.ZH.X.A.XT*  AT. Z) 

042  MT«NP>2 

043  N*24MT4| 

044  MR  I TE 1 3. 29) IZI J ) . 3»l *N) 

045  29  FORMAT I  *  2* /I IX. 6E1 1 .41* 

046  CALL  PLANE 1 1.1*1 *NP.NT .RH. ZH.XT. AT *THR*R| 

047  00  22  J* 1 .NT 

046  B!J)«R!J) 

049  JI*J«NT 

050  61 Jl )*-R IJI ) 

051  22  CONTINUE 

052  6INI«-AINI 

053  WRITE I3.231I8IJ) » J*l *N) 

054  23  FORMAT! •  B* /I I X.6E 1 I .4) ) 

055  CALL  OECOMPIN.IPS.Z) 

056  CALL  SOLVE! N*IPS.Z*B*C» 

057  U*<0..1.) 

058  MRITEI3.21) 

059  21  FORMAT!  •  REAL  JT  (MAG  JT  MAG  JT*» 

060  00  24  Jal.MT 
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061  CI*2./RH(J»||(C(J) 

062  C2»CA0S(CI) 

063  9RITE13.2S1  CI.C2 

064  25  FORMAT IIX.3Ett.4l 

065  24  CONTINUE 

066  NRITE13.26> 

067  26  FORMAT!  •  REA*.  JP  I  MAC  JR  MAS  JR*  | 

066  NP«NR>I 

069  00  27  J>I,NP 

070  CI*4./<RH<J»4RH< J*ll )«U4C( J4MT1 

071  C2sCABS(CII 

072  WR1TEI3.251  C1.C2 

073  27  CONTINUE 

074  STOP 

075  ENO 

•  DATA 
2  20 

-  0.5 7735 03E+00  0.5773S03E400 
0  *  t  900000E+01  0.  10000006401 

-  0. 993 1 2866400-0.963971 96400-0 .9 1223446400-0 *6391 1706400-0 .74033196*00 

-  0.6360537E400- 0.5 1066 7064 00-0. 3 73 7061 6400-0. 22778S9E400-0 .76526 526-01 
0. 76S2652E— 01  0.22778596*00  0.37370616400  0.51086706*00  0.63605376*00 
0.74633196*00  0.6391 1  TOE *- 00  0.9 122344E400  0.963971 9E*00  0.99312866400 
0. 176 I40IE-0I  0 .4060 1 4  3E— 0  I  0.62672056-01  0. 6327674E- 01  O • 10 193016400 
0.  1  18  1945E.OO  0.  13166866400  0.142096164-00  0 . I 491 730E* 00  0 . 1S27534E4-00 
0.  15275346*00  0. 1491  730E4-00  0.142096164-00  0.  1 3 I6886E* 00  0.116 19456*00 
0*101  9301  E-400  C.8327674E-01  0  .62672056-01  0.4  060143E-01  O  .  1 76 1 401C-01 
11  0.  12566376*01  0. 0000000E4-00 

0.0000  0.2334  0.4540  0.6494  0.6090  0.9239  0.9877  0.9969  0.9511  0.6526 

0.7071 

-1.0000  -0.9724  -0.8910  -0.7604  -0.5878  -0.3827  -0.1564  0.0765  0.3090  0.5225 

0.7071 
4ST0P 
/* 

// 

PRINTED  OUTPUT 
NT  NPMI 
2  20 
XT 

-0.57735036*00  O.S773503E400 
AT 

0. 1000000E401  0.1000009E4-01 
X 

-  0.993 12866 4-00—  0.963971 9E4 00— 0  .9 1 22344 E+ 00— 0.0391 I 70E*00— 0.74633196*00 

-  0.63605376*00— 0.51086706*00— 0.3737061 6* 00-0.22776596* 00—0. 7652652E— 01 
O. 76526526—01  0.22770S9E400  0.3737061  E 4- 00  0.5I00670E*00  0.636053764-00 
0.74633196400  0.6391 1 70E4  00  0.91223446400  0.96397I9E400  0.99312866400 

A 

0.1 76  1 4  01 E— 01  0.4060 1 43E— 01  0.62672066-01  0.63276756-01  0.10193016400 
0.1 I8194SE400  0. 131 66866400  0 • 1420961E400  0.14917306400  0. 15275346400 
0. 15275346400  0.14917306400  0. 1420901 6400  0 .13160066400  0 .1 1819456400 
0.10193016400  0.O327O7SE— 01  0.62672086-01  0.4060143E-01  0.17614016-01 
NP  BK  THR 

II  0.12560376401  0.00000006400 


RH 

0.0000 

0.7071 

ZH 

0.2334  0.4540 

0.6494 

0.8090 

0.0239 

0.9077 

0.9909 

0.9511 

0.0520 

-I. 0000 
0.7071 

-0.9724  -0.0910 

-0.7604 

-0.  5870 

-0.3027 

-0.1504 

0.0705 

0.3090 

0.5225 

0.83636-01-0. 97786401  0.7S52E-01  0.29746401  t.«]*4E-OI  0*77236400 
0*49196— 01  0* 28936400  0*34806-01  0.14916400  0.21266-01  0. 91276-01 
0*10246-01  0*64776—01  0*19546-02  0.5293C-0I-0.3673E-02  0.48426-01 
0.13486402  0.67016-01  0*24136401  0*84146-01  0*42506400  0*78586-01 
0.13696400  0.71036-01  0*76336-01  0*62136-01  0.S943E-0I  0*52656-01 
0*54376—01  0*43266-01  0*52756-01  0.34SSE-01  6*52826-01  6*26846-61 
0*51346—01  0*20636-01 

0*31516400-0.83856400  0*36316400-0.73526400  0*40466406-0*56666468 
0. 39666400-0. 36S6E400  0.30746400-0*16696400  0* I384E40O-0.37866-0I 

-  0  *  70376-01-0. 1 74 86- 0 1-0*  25876 4 00- 0* 1 1 52 6400-0*3766640 0-0*  28636406 
0*87556400  0.30706400  0*85276406  0*36606400  0*79596400  0*47496466 
0*69176400  0.61616400  0*52796400  0.76026400  0*30606400  0.6730E400 
0*45036-01  0*«236E400-0*22196400  0.69756400-0*45966400  0.00316466 

-0*64386400  0.66536400 

REAL  JT  I MAO  J7  MAC  JT 

0.11426401  0.11966401  0*16556401 
0.81876400  0.11996401  0.14SIE401 
0*34326400  0.1 163640 1  0* (2126401 
-0.2206E400  0.10356401  0*10586401 
-0.77666400  0.78276400  0*11036401 
-0*12226401  0*41616400  0.12916401 
-0*14746401-0.61466-02  0*14746401 
-0.14896401-0*38776400  0*15386401 
-0.12486401-0.61336400  0*13916401 
REAL  JP  (MAC  JP  MAC  JP 

-0.12096401-0. 11846401  0*16926401 
-0*10766401-0.10946401  0*15346401 
-0.94866400-0.94186400  0*13376401 
-0.85846400-0.72136400  0.11216401 

-  0.8928E4Q 0—0 *48636400  0*10176401 
-0*11026401-0*33416400  0*11526401 

-  0*14806401-0. 3804E400  0.15286401 
-0.20026401-0.73556400  0*21336401 
-0*26296401-0*14336401  0*29956401 
-0.63816401-0.51796401  0*82196401 
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